Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jakub.both/dune-istl
  • eduardo.bueno/dune-istl
  • tkoch/dune-istl
  • pipping/dune-istl
  • andreas.brostrom/dune-istl
  • stephan.hilb/dune-istl
  • max.kahnt/dune-istl
  • patrick.jaap/dune-istl
  • andreas.thune/dune-istl
  • dominic/dune-istl
  • ansgar/dune-istl
  • exadune/dune-istl
  • lars.lubkoll/dune-istl
  • govind.sahai/dune-istl
  • maikel.nadolski/dune-istl
  • markus.blatt/dune-istl
  • core/dune-istl
  • lisa_julia.nebel/dune-istl
  • michael.sghaier/dune-istl
  • liesel.schumacher/dune-istl
  • Xinyun.Li/dune-istl
  • lukas.renelt/dune-istl
  • simon.praetorius/dune-istl
  • rene.milk/dune-istl
  • lasse.hinrichsen/dune-istl
  • nils.dreier/dune-istl
  • claus-justus.heine/dune-istl
  • felix.mueller/dune-istl
  • kilian.weishaupt/dune-istl
  • lorenzo.cerrone/dune-istl
  • jakob.torben/dune-istl
  • liam.keegan/dune-istl
  • alexander.mueller/dune-istl
33 results
Show changes
Commits on Source (215)
Showing
with 611 additions and 243 deletions
......@@ -22,33 +22,6 @@ variables:
# than one thread.
OMP_NUM_THREADS: 1
debian-11-gcc-9-17-with-checking:
image: registry.dune-project.org/docker/ci/debian:11
script: duneci-standard-test
variables:
DUNECI_TOOLCHAIN: gcc-9-17
DUNECI_CMAKE_FLAGS: "CC=gcc-9 CXX=g++-9 -DCMAKE_CXX_FLAGS='-std=c++17 -O2 -g -Wall -fdiagnostics-color=always -DDUNE_ISTL_WITH_CHECKING' -DDUNE_MAX_TEST_CORES=4 -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
tags: [duneci]
# no numpy version available for this image on the gitlab registry
# disable python bindings for now
ubuntu:18.04 clang-6-17:
image: registry.dune-project.org/docker/ci/ubuntu:18.04
script: duneci-standard-test
stage: test
variables:
DUNECI_TOOLCHAIN: clang-6-17
DUNECI_CMAKE_FLAGS: "-DDUNE_ENABLE_PYTHONBINDINGS=OFF"
tags: [duneci]
ubuntu:18.04 clang-5-17:
image: registry.dune-project.org/docker/ci/ubuntu:18.04
script: duneci-standard-test
stage: test
variables:
DUNECI_TOOLCHAIN: clang-5-17
DUNECI_CMAKE_FLAGS: "-DDUNE_ENABLE_PYTHONBINDINGS=OFF"
tags: [duneci]
# Check for spelling mistakes in text
code-spelling-check:
stage: .pre
......
......@@ -3,8 +3,43 @@ SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE
SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
-->
# Master (will become release 2.10)
# Master (will become release 2.11)
- `FastAMG` can now be used with non-nested matrices of type `BCRSMatrix<T>`
where `T` is a plain number type instead of `FieldMatrix<T,m,n>`.
# Release 2.10
- Improve testing support on Laplacian matrices with an optional diagonal regularization parameter.
- Base the implementation of `VariableBlockVector` on `std::vector` as the storage type. Note that
this prevents from using `bool` as block type that was possible before.
- A method `BCRSMatrix::setIndicesNoSort()` was added. Similar
to `BCRSMatrix::setIndices()` this allows to insert all indices
of a row at once, but - incontrast to the latter - does not sort them.
This allows to speed up insertion if indices are already sorted.
- `UMFPACK` is extended to arbitrary blocked matrix structures. This includes `MultiTypeBlockMatrix`. The external interface is unchanged.
However, internally the `BCCSMatrixInitializer` is replaced by direct calls of `flatMatrixForEach` similar to `Cholmod`. This requires a
compatible vector field of "ignored" degrees of freedom. The method `setSubMatrix` with a top-level index set is preserved for compatibility.
- The internal storage in `MatrixIndexSet` was changed from `std::set` to a sorted `std::vector`
(with a `std::set` fallback for very dense rows) to improve performance. The stored index
type was changed from `std::size_t` to `uint32_t` to reduce memory consumption and improve
performance. Hence, `MatrixIndexSet` can no longer be used for very large matrices with more
than 2^32 columns.
- Added flag 'useFixedOrder' to the coarsen method of AMGs ParallelIndicesCoarsener.
If set to true, during the creation of the coarser matrix (by accumulation and restriction
to fewer processes) the coarse matrix rows are ordered in a fixed manner, making parallel runs
reproducible; the runtime is possibly not ideal though.
If set to false (which is the default), the row order depends on the order of messages received from the
processes responsible for the respective parts of the finer grid. Then the indices on the coarser grid
may differ from run to run.
- Define `field_type` and `real_type` in `MultiTypeBlock[Matrix|Vector]` only if a common type of these
types exist over all blocks in the container. Otherwise it is defined as `Std::nonesuch`.
## Deprecations and removals
......@@ -21,6 +56,9 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
- Removed deprecated `MultyTypeBlockVector::count()`, use `N()` instead.
- Deprecated `writeSVGMatrix` with `std::stream` as the second argument. Use the method
with `std::stream` as the first argument.
# Release 2.9
......@@ -48,6 +86,11 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
used by `CHOLMOD` itself to store the Cholesky factor. This can be used to
use the more advanced features of `CHOLMOD`.
- The `Cholmod` class now has an additional template parameter `Index`, which specifies the
type used by `CHOLMOD` for integers. The types `int` (the default) and `SuiteSparse_long`
(which is usually `long int`) are supported. Use `SuiteSparse_long` if your problems
are too big for `int`.
- You can now multiply objects of type `ScaledIdentityMatrix` by scalars
using `operator*`.
......
# SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
# set up project
cmake_minimum_required(VERSION 3.16)
project(dune-istl C CXX)
# general stuff
cmake_minimum_required(VERSION 3.13)
# guess build tree of dune-common
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
......@@ -23,9 +20,20 @@ list(APPEND CMAKE_MODULE_PATH ${dune-common_MODULE_PATH}
#include the dune macros
include(DuneMacros)
# deactivate global include-directories
dune_policy(SET DP_DEFAULT_INCLUDE_DIRS NEW)
# start a dune project with information from dune.module
dune_project()
# create library target and export it as Dune::ISTL
dune_add_library(duneistl INTERFACE
EXPORT_NAME ISTL
LINK_LIBRARIES Dune::Common)
# set include directories to target
dune_default_include_directories(duneistl INTERFACE)
add_subdirectory(cmake/modules)
add_subdirectory(dune)
add_subdirectory(doc)
......@@ -37,4 +45,4 @@ if( DUNE_ENABLE_PYTHONBINDINGS )
endif()
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
finalize_dune_project()
......@@ -3,37 +3,41 @@ Copyright holders:
2014--2016 Marco Agnese
2003--2013 Peter Bastian
2020 Jean Benezech
2004--2022 Markus Blatt
2014--2022 Ansgar Burchardt
2004--2024 Markus Blatt
2022 Eduardo Bueno
2022 Samuel Burbulla
2014--2023 Ansgar Burchardt
2004--2005 Adrian Burri
2017--2018 Matthew Collins
2008--2021 Andreas Dedner
2008--2024 Andreas Dedner
2018--2022 Nils-Arne Dreier
2003--2022 Christian Engwer
2003--2023 Christian Engwer
2005--2019 Jorrit Fahlke
2008--2017 Bernd Flemisch
2017--2020 Janick Gerstenberger
2015--2020 Felix Gruber
2005--2021 Carsten Gräser
2012--2021 Christoph Grüninger
2005--2024 Carsten Gräser
2012--2024 Christoph Grüninger
2016--2020 René Heß
2018 Claus-Justus Heine
2016--2019 Stephan Hilb
2018 Lasse Hinrichsen
2018--2024 Lasse Hinrichsen
2012--2013 Olaf Ippisch
2018--2021 Patrick Jaap
2018--2024 Patrick Jaap
2013--2019 Dominic Kempf
2015 Emmanouil Kiagias
2004--2022 Robert Klöfkorn
2016--2022 Timo Koch
2004--2023 Robert Klöfkorn
2016--2023 Timo Koch
2007 Sreejith Pulloor Kuttanikkad
2013--2016 Arne Morten Kvarving
2009--2013 Andreas Lauser
2012--2017 Tobias Malkmus
2007--2009 Sven Marnach
2013 René Milk
2024 Alex Müller
2013--2019 Steffen Müthing
2016 Maikel Nadolski
2023--2024 Lisa Julia Nebel
2003--2005 Thimo Neubauer
2010--2012 Rebecca Neumann
2008--2018 Martin Nolte
......@@ -41,15 +45,17 @@ Copyright holders:
2013--2014 Marian Piatkowski
2011--2016 Elias Pipping
2013 Jurgis Pods
2018--2022 Simon Praetorius
2018--2024 Simon Praetorius
2009 Atgeirr Rasmussen
2018--2022 Lukas Renelt
2021--2022 Santiago Ospina De Los Ríos
2021--2024 Santiago Ospina De Los Ríos
2008--2013 Uli Sack
2004--2021 Oliver Sander
2023 Vaishnavi Sanchi
2004--2024 Oliver Sander
2015--2019 Linus Seelinger
2013 Bård Skaflestad
2018 Mathis Springwald
2023 Jakob Torben
2006--2010 Martin Weiser
2019--2020 Kilian Weishaupt
2015 Sebastian Westerheide
......
......@@ -12,12 +12,19 @@
#
# A list of targets to use ARPACKPP with.
#
include_guard(GLOBAL)
set_package_properties("ARPACK" PROPERTIES
PURPOSE "Solve large scale eigenvalue problems")
set_package_properties("ARPACKPP" PROPERTIES
PURPOSE "C++ interface for ARPACK")
function(add_dune_arpackpp_flags _targets)
if(ARPACKPP_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} PUBLIC ${ARPACKPP_DUNE_LIBRARIES})
target_compile_definitions(${_target} PUBLIC ENABLE_ARPACKPP=1)
target_compile_definitions(${_target} PUBLIC HAVE_ARPACKPP=1)
target_compile_options(${_target} PUBLIC ${ARPACKPP_DUNE_COMPILE_FLAGS})
endforeach()
endif()
......
......@@ -12,6 +12,10 @@
#
# A list of targets to use SuperLU with.
#
include_guard(GLOBAL)
set_package_properties("SuperLU" PROPERTIES
PURPOSE "Direct solver for linear system, based on LU decomposition")
# set HAVE_SUPERLU for config.h
set(HAVE_SUPERLU ${SuperLU_FOUND})
......@@ -19,7 +23,7 @@ set(HAVE_SUPERLU ${SuperLU_FOUND})
# register all SuperLU related flags
if(SuperLU_FOUND)
dune_register_package_flags(
COMPILE_DEFINITIONS "ENABLE_SUPERLU=1"
COMPILE_DEFINITIONS "HAVE_SUPERLU=1"
LIBRARIES SuperLU::SuperLU)
endif()
......@@ -28,7 +32,7 @@ function(add_dune_superlu_flags _targets)
if(SuperLU_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} PUBLIC SuperLU::SuperLU)
target_compile_definitions(${_target} PUBLIC ENABLE_SUPERLU=1)
target_compile_definitions(${_target} PUBLIC HAVE_SUPERLU=1)
endforeach()
endif()
endfunction(add_dune_superlu_flags)
......
......@@ -13,8 +13,6 @@ find_package(SuperLU 5.0)
include(AddSuperLUFlags)
find_package(ARPACKPP)
include(AddARPACKPPFlags)
find_package(SuiteSparse OPTIONAL_COMPONENTS CHOLMOD LDL SPQR UMFPACK)
include(AddSuiteSparseFlags)
# enable / disable backwards compatibility w.r.t. category
set(DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE 1
......
......@@ -25,6 +25,11 @@
# system paths.
#
# text for feature summary
set_package_properties("ARPACK" PROPERTIES
URL "https://www.arpack.org"
DESCRIPTION "ARnoldi PACKage")
# look for library, only at positions given by the user
find_library(ARPACK_LIBRARY
NAMES "arpack"
......@@ -81,8 +86,3 @@ else()
"Determining location of ARPACK failed:\n"
"Libraries to link against: ${ARPACK_LIBRARIES}\n\n")
endif()
# text for feature summary
set_package_properties("ARPACK" PROPERTIES
DESCRIPTION "ARnoldi PACKage"
PURPOSE "Solve large scale eigenvalue problems")
......@@ -28,6 +28,11 @@
# system paths.
#
# text for feature summary
set_package_properties("ARPACKPP" PROPERTIES
URL "https://github.com/m-reuter/arpackpp"
DESCRIPTION "ARPACK++")
# find ARPACK which is required by ARPACK++
find_package(ARPACK)
......@@ -123,17 +128,12 @@ else()
"Libraries to link against: ${ARPACKPP_LIBRARIES}\n\n")
endif()
# set HAVE_ARPACKPP for config.h
# set HAVE_ARPACKPP for conditional tests
set(HAVE_ARPACKPP ${ARPACKPP_FOUND})
# register all ARPACK++ related flags
if(ARPACKPP_FOUND)
dune_register_package_flags(COMPILE_DEFINITIONS "ENABLE_ARPACKPP=1"
dune_register_package_flags(COMPILE_DEFINITIONS "HAVE_ARPACKPP=1"
LIBRARIES "${ARPACKPP_LIBRARIES}"
COMPILE_OPTIONS "${ARPACKPP_DUNE_COMPILE_FLAGS}")
endif()
# text for feature summary
set_package_properties("ARPACKPP" PROPERTIES
DESCRIPTION "ARPACK++"
PURPOSE "C++ interface for ARPACK")
......@@ -25,8 +25,8 @@
# text for feature summary
include(FeatureSummary)
set_package_properties("SuperLU" PROPERTIES
DESCRIPTION "Supernodal LU"
PURPOSE "Direct solver for linear system, based on LU decomposition")
URL "https://portal.nersc.gov/project/sparse/superlu/"
DESCRIPTION "Supernodal LU")
set(SUPERLU_INT_TYPE "int" CACHE STRING
"The integer version that SuperLU was compiled for (Default is int.
......
......@@ -30,16 +30,10 @@
/* end private */
/* Define to ENABLE_SUPERLU if the SuperLU library is available */
#cmakedefine HAVE_SUPERLU ENABLE_SUPERLU
/* Define to the integer type that SuperLU was compiled for
See e.g. what int_t is defined to in slu_sdefs.h */
#cmakedefine SUPERLU_INT_TYPE @SUPERLU_INT_TYPE@
/* Define to ENABLE_ARPACKPP if the ARPACK++ library is available */
#cmakedefine HAVE_ARPACKPP ENABLE_ARPACKPP
/* Define to the version of dune-istl */
#define DUNE_ISTL_VERSION "${DUNE_ISTL_VERSION}"
......
......@@ -2,11 +2,11 @@
# SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
Module: dune-istl
Version: 2.10-git
Version: 2.11-git
Author: The Dune Core developers
Maintainer: dune-devel@lists.dune-project.org
Description: Iterative solver template library which provides generic sparse matrix/vector classes and a variety of solvers based on these classes
URL: https://gitlab.dune-project.org/core/dune-istl
Python-Requires:
Depends: dune-common (>= 2.10)
Depends: dune-common (>= 2.11)
Whitespace-Hook: Yes
......@@ -18,6 +18,7 @@ install(FILES
btdmatrix.hh
bvector.hh
cholmod.hh
dilu.hh
foreach.hh
gsetc.hh
ildl.hh
......
......@@ -5,11 +5,12 @@
#ifndef DUNE_ISTL_BASEARRAY_HH
#define DUNE_ISTL_BASEARRAY_HH
#include "assert.h"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <memory>
#include <algorithm>
#include <type_traits>
#include "istlexception.hh"
#include <dune/common/iteratorfacades.hh>
......@@ -41,13 +42,9 @@ namespace Imp {
Setting the compile time switch DUNE_ISTL_WITH_CHECKING
enables error checking.
\todo There shouldn't be an allocator argument here, because the array is 'unmanaged'.
And indeed, of the allocator, only its size_type is used. Hence, the signature
of this class should be changed to <class B, int stype>
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
template<class B, class ST=std::size_t >
class base_array_unmanaged
{
public:
......@@ -57,11 +54,8 @@ namespace Imp {
//! export the type representing the components
typedef B member_type;
//! export the allocator type
typedef A allocator_type;
//! the type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! the type used for references
using reference = B&;
......@@ -104,17 +98,27 @@ namespace Imp {
friend class RealIterator<ValueType>;
//! constructor
RealIterator ()
: p(0), i(0)
{}
RealIterator () = default;
RealIterator (const B* _p, B* _i) : p(_p), i(_i)
{ }
RealIterator (const B* _p, B* _i)
: p(_p), i(_i)
{}
RealIterator(const RealIterator<ValueType>& it)
: p(it.p), i(it.i)
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator (const RealIterator<T_>& other)
: p(other.p), i(other.i)
{}
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator& operator= (const RealIterator<T_>& other)
{
p = other.p;
i = other.i;
return *this;
}
//! return index
size_type index () const
{
......@@ -170,8 +174,8 @@ namespace Imp {
i+=d;
}
const B* p;
B* i;
const B* p = nullptr;
B* i = nullptr;
};
//! iterator type for sequential access
......@@ -302,7 +306,7 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
template<class B, class ST=std::size_t >
class compressed_base_array_unmanaged
{
public:
......@@ -312,11 +316,8 @@ namespace Imp {
//! export the type representing the components
typedef B member_type;
//! export the allocator type
typedef A allocator_type;
//! The type used for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! the type used for references
using reference = B&;
......@@ -359,22 +360,29 @@ namespace Imp {
friend class RealIterator<ValueType>;
//! constructor
RealIterator ()
: p(0), j(0), i(0)
{}
RealIterator () = default;
//! constructor
RealIterator (B* _p, size_type* _j, size_type _i)
: p(_p), j(_j), i(_i)
{ }
{}
/**
* @brief Copy constructor from mutable iterator
*/
RealIterator(const RealIterator<ValueType>& it)
: p(it.p), j(it.j), i(it.i)
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator (const RealIterator<T_>& other)
: p(other.p), j(other.j), i(other.i)
{}
template <class T_,
std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
RealIterator& operator= (const RealIterator<T_>& other)
{
p = other.p;
j = other.j;
i = other.i;
return *this;
}
//! equality
bool equals (const RealIterator<ValueType>& it) const
......@@ -434,9 +442,9 @@ namespace Imp {
return p[i];
}
B* p;
size_type* j;
size_type i;
B* p = nullptr;
size_type* j = nullptr;
size_type i = 0;
};
/** @brief The iterator type. */
......
......@@ -493,12 +493,12 @@ namespace Dune {
//! export the allocator type
typedef A allocator_type;
//! implement row_type with compressed vector
typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
//! The type for the index access and the size
typedef typename A::size_type size_type;
//! implement row_type with compressed vector
typedef Imp::CompressedBlockVectorWindow<B,size_type> row_type;
//! The type for the statistics object returned by compress()
typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
......@@ -775,6 +775,7 @@ namespace Dune {
* @param _avg expected average number of entries per row
* @param compressionBufferSize fraction of _n*_avg which is expected to be
* needed for elements that exceed _avg entries per row.
* \param bm Matrix build mode
*
*/
BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
......@@ -1218,11 +1219,35 @@ namespace Dune {
*pos = col;
}
//! Set all column indices for row from the given iterator range.
/**
* The iterator range has to be of the same length as the previously set row size.
* The entries in the iterator range must be sorted and must not contain duplicate values.
* This method will insert the indices in the given order.
*
* Calling this method overwrites any previously set column indices!
*
*/
template<typename It>
void setIndicesNoSort(size_type row, It begin, It end)
{
size_type row_size = r[row].size();
size_type* col_begin = r[row].getindexptr();
size_type* col_end;
// consistency check between allocated row size and number of passed column indices
if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
DUNE_THROW(BCRSMatrixError,"Given size of row " << row
<< " (" << row_size
<< ") does not match number of passed entries (" << (col_end - col_begin) << ")");
}
//! Set all column indices for row from the given iterator range.
/**
* The iterator range has to be of the same length as the previously set row size.
* The entries in the iterator range do not have to be in any particular order, but
* must not contain duplicate values.
* This method will insert the indices and sort them afterwards.
*
* Calling this method overwrites any previously set column indices!
*/
......@@ -2199,12 +2224,13 @@ namespace Dune {
* setup matrix) results in all the structure and data being lost. Please
* call deallocate() before calling allocate in this case.
*
* @param row The number of rows the matrix should contain.
* @param rows The number of rows the matrix should contain.
* @param columns the number of columns the matrix should contain.
* @param allocationSize The number of nonzero entries the matrix should hold (if omitted
* defaults to 0).
* @param allocateRow Whether we have to allocate the row pointers, too. (Defaults to
* @param allocateRows Whether we have to allocate the row pointers, too. (Defaults to
* true)
* \param allocate_data Whether data array needs to be allocated
*/
void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
{
......
......@@ -80,8 +80,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class block_vector_unmanaged : public base_array_unmanaged<B,A>
template<class B, class ST=std::size_t >
class block_vector_unmanaged : public base_array_unmanaged<B,ST>
{
public:
......@@ -91,17 +91,14 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! The size type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! make iterators available as types
typedef typename base_array_unmanaged<B,A>::iterator Iterator;
typedef typename base_array_unmanaged<B,ST>::iterator Iterator;
//! make iterators available as types
typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
typedef typename base_array_unmanaged<B,ST>::const_iterator ConstIterator;
//! for STL compatibility
typedef B value_type;
......@@ -177,8 +174,8 @@ namespace Imp {
* @param y other (compatible) vector
* @return
*/
template<class OtherB, class OtherA>
auto operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
template<class OtherB, class OtherST>
auto operator* (const block_vector_unmanaged<OtherB,OtherST>& y) const
{
typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
PromotedType sum(0);
......@@ -198,8 +195,8 @@ namespace Imp {
* @param y other (compatible) vector
* @return
*/
template<class OtherB, class OtherA>
auto dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
template<class OtherB, class OtherST>
auto dot(const block_vector_unmanaged<OtherB,OtherST>& y) const
{
typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
PromotedType sum(0);
......@@ -338,7 +335,7 @@ namespace Imp {
protected:
//! make constructor protected, so only derived classes can be instantiated
block_vector_unmanaged () : base_array_unmanaged<B,A>()
block_vector_unmanaged () : base_array_unmanaged<B,ST>()
{ }
};
......@@ -391,7 +388,7 @@ namespace Imp {
enables error checking.
*/
template<class B, class A=std::allocator<B> >
class BlockVector : public Imp::block_vector_unmanaged<B,A>
class BlockVector : public Imp::block_vector_unmanaged<B,typename A::size_type>
{
public:
......@@ -410,10 +407,10 @@ namespace Imp {
typedef typename A::size_type size_type;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
//===== constructors and such
......@@ -552,7 +549,7 @@ namespace Imp {
BlockVector& operator= (const field_type& k)
{
// forward to operator= in base class
(static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
return *this;
}
......@@ -619,7 +616,7 @@ namespace Imp {
#else
template<class B, class A=std::allocator<B> >
#endif
class BlockVectorWindow : public Imp::block_vector_unmanaged<B,A>
class BlockVectorWindow : public Imp::block_vector_unmanaged<B,typename A::size_type>
{
public:
......@@ -638,15 +635,15 @@ namespace Imp {
typedef typename A::size_type size_type;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
//! make iterators available as types
typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
//===== constructors and such
//! makes empty array
BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
BlockVectorWindow () : Imp::block_vector_unmanaged<B,size_type>()
{ }
//! make array from given pointer and size
......@@ -682,7 +679,7 @@ namespace Imp {
//! assign from scalar
BlockVectorWindow& operator= (const field_type& k)
{
(static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
return *this;
}
......@@ -741,8 +738,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,A>
template<class B, class ST=std::size_t >
class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,ST>
{
public:
......@@ -754,17 +751,14 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! make iterators available as types
typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
typedef typename compressed_base_array_unmanaged<B,ST>::iterator Iterator;
//! make iterators available as types
typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
typedef typename compressed_base_array_unmanaged<B,ST>::const_iterator ConstIterator;
//! The type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//===== assignment from scalar
......@@ -961,7 +955,7 @@ namespace Imp {
protected:
//! make constructor protected, so only derived classes can be instantiated
compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,ST>()
{ }
//! return true if index sets coincide
......@@ -995,8 +989,8 @@ namespace Imp {
\internal This class is an implementation detail, and should not be used outside of dune-istl.
*/
template<class B, class A=std::allocator<B> >
class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,A>
template<class B, class ST=std::size_t >
class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,ST>
{
public:
......@@ -1008,22 +1002,19 @@ namespace Imp {
//! export the type representing the components
typedef B block_type;
//! export the allocator type
typedef A allocator_type;
//! The type for the index access
typedef typename A::size_type size_type;
typedef ST size_type;
//! make iterators available as types
typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
typedef typename compressed_block_vector_unmanaged<B,ST>::Iterator Iterator;
//! make iterators available as types
typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
typedef typename compressed_block_vector_unmanaged<B,ST>::ConstIterator ConstIterator;
//===== constructors and such
//! makes empty array
CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,ST>()
{ }
//! make array from given pointers and size
......@@ -1062,7 +1053,7 @@ namespace Imp {
//! assign from scalar
CompressedBlockVectorWindow& operator= (const field_type& k)
{
(static_cast<compressed_block_vector_unmanaged<B,A>&>(*this)) = k;
(static_cast<compressed_block_vector_unmanaged<B,ST>&>(*this)) = k;
return *this;
}
......
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#pragma once
#if HAVE_SUITESPARSE_CHOLMOD
#if HAVE_SUITESPARSE_CHOLMOD || defined DOXYGEN
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include<dune/istl/solver.hh>
#include <dune/istl/solverfactory.hh>
#include <dune/istl/solverregistry.hh>
#include <dune/istl/foreach.hh>
#include <vector>
......@@ -66,16 +68,157 @@ namespace Impl{
});
}
// wrapper class for C function calls to CHOLMOD itself.
// The CHOLMOD API has different functions for different index types.
template <class Index>
struct CholmodMethodChooser;
// specialization using 'int' to store indices
template <>
struct CholmodMethodChooser<int>
{
[[nodiscard]]
static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
{
return ::cholmod_allocate_dense(nrow,ncol,d,xtype,c);
}
[[nodiscard]]
static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
{
return ::cholmod_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
}
[[nodiscard]]
static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
{
return ::cholmod_analyze(A,c);
}
static int defaults(cholmod_common *c)
{
return ::cholmod_defaults(c);
}
static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
{
return ::cholmod_factorize(A,L,c);
}
static int finish(cholmod_common *c)
{
return ::cholmod_finish(c);
}
static int free_dense (cholmod_dense **X, cholmod_common *c)
{
return ::cholmod_free_dense(X,c);
}
static int free_factor(cholmod_factor **L, cholmod_common *c)
{
return ::cholmod_free_factor(L,c);
}
static int free_sparse(cholmod_sparse **A, cholmod_common *c)
{
return ::cholmod_free_sparse(A,c);
}
[[nodiscard]]
static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
{
return ::cholmod_solve(sys,L,B,c);
}
static int start(cholmod_common *c)
{
return ::cholmod_start(c);
}
};
// specialization using 'SuiteSparse_long' to store indices
template <>
struct CholmodMethodChooser<SuiteSparse_long>
{
[[nodiscard]]
static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
{
return ::cholmod_l_allocate_dense(nrow,ncol,d,xtype,c);
}
[[nodiscard]]
static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
{
return ::cholmod_l_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
}
[[nodiscard]]
static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
{
return ::cholmod_l_analyze(A,c);
}
static int defaults(cholmod_common *c)
{
return ::cholmod_l_defaults(c);
}
static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
{
return ::cholmod_l_factorize(A,L,c);
}
static int finish(cholmod_common *c)
{
return ::cholmod_l_finish(c);
}
static int free_dense (cholmod_dense **X, cholmod_common *c)
{
return ::cholmod_l_free_dense(X,c);
}
static int free_factor (cholmod_factor **L, cholmod_common *c)
{
return ::cholmod_l_free_factor(L,c);
}
static int free_sparse(cholmod_sparse **A, cholmod_common *c)
{
return ::cholmod_l_free_sparse(A,c);
}
[[nodiscard]]
static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
{
return ::cholmod_l_solve(sys,L,B,c);
}
static int start(cholmod_common *c)
{
return ::cholmod_l_start(c);
}
};
} //namespace Impl
/** @brief Dune wrapper for SuiteSparse/CHOLMOD solver
*
* This class implements an InverseOperator between Vector types
*
* \tparam Vector Data type for solution and right-hand-side vectors
* \tparam Index Type used by CHOLMOD for indices. Must be either 'int' or 'SuiteSparse_long'
* (which is usually `long int`).
*/
template<class Vector>
template<class Vector, class Index=int>
class Cholmod : public InverseOperator<Vector, Vector>
{
static_assert(std::is_same_v<Index,int> || std::is_same_v<Index,SuiteSparse_long>,
"Index type must be either 'int' or 'SuiteSparse_long'!");
using CholmodMethod = Impl::CholmodMethodChooser<Index>;
public:
/** @brief Default constructor
......@@ -85,7 +228,7 @@ public:
*/
Cholmod()
{
cholmod_start(&c_);
CholmodMethod::start(&c_);
}
/** @brief Destructor
......@@ -96,8 +239,8 @@ public:
~Cholmod()
{
if (L_)
cholmod_free_factor(&L_, &c_);
cholmod_finish(&c_);
CholmodMethod::free_factor(&L_, &c_);
CholmodMethod::finish(&c_);
}
// forbid copying to avoid freeing memory twice
......@@ -107,7 +250,7 @@ public:
/** @brief simple forward to apply(X&, Y&, InverseOperatorResult&)
*/
void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
{
apply(x,b,res);
}
......@@ -117,7 +260,7 @@ public:
* The method assumes that setMatrix() was called before
* In the case of a given ignore field the corresponding entries of both in x and b will stay untouched in this method.
*/
void apply(Vector& x, Vector& b, InverseOperatorResult& res)
void apply(Vector& x, Vector& b, InverseOperatorResult& res) override
{
// do nothing if N=0
if ( nIsZero_ )
......@@ -144,13 +287,14 @@ public:
});
// create a cholmod dense object
auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
auto b3 = make_cholmod_dense(CholmodMethod::allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
// cast because void-ptr
auto b4 = static_cast<double*>(b3->x);
std::copy(b2.get(), b2.get() + L_->n, b4);
// solve for a cholmod x object
auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
auto x3 = make_cholmod_dense(CholmodMethod::solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
// cast because void-ptr
auto xp = static_cast<double*>(x3->x);
......@@ -199,8 +343,8 @@ public:
void setMatrix(const Matrix& matrix, const Ignore* ignore)
{
// count the number of entries and diagonal entries
int nonZeros = 0;
int numberOfIgnoredDofs = 0;
size_t nonZeros = 0;
size_t numberOfIgnoredDofs = 0;
auto [flatRows,flatCols] = flatMatrixForEach( matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
......@@ -216,38 +360,38 @@ public:
numberOfIgnoredDofs = std::count(flatIgnore.begin(),flatIgnore.end(),true);
}
// Total number of rows
int N = flatRows - numberOfIgnoredDofs;
nIsZero_ = (N <= 0);
nIsZero_ = (size_t(flatRows) <= numberOfIgnoredDofs);
if ( nIsZero_ )
{
return;
}
// Total number of rows
size_t N = flatRows - numberOfIgnoredDofs;
/*
* CHOLMOD uses compressed-column sparse matrices, but for symmetric
* matrices this is the same as the compressed-row sparse matrix used
* by DUNE. So we can just store Mᵀ instead of M (as M = Mᵀ).
*/
const auto deleter = [c = &this->c_](auto* p) {
cholmod_free_sparse(&p, c);
CholmodMethod::free_sparse(&p, c);
};
auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
cholmod_allocate_sparse(N, // # rows
N, // # cols
nonZeros, // # of nonzeroes
1, // indices are sorted ( 1 = true)
1, // matrix is "packed" ( 1 = true)
-1, // stype of matrix ( -1 = consider the lower part only )
CHOLMOD_REAL, // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
&c_ // cholmod_common ptr
), deleter);
CholmodMethod::allocate_sparse(N, // # rows
N, // # cols
nonZeros, // # of nonzeroes
1, // indices are sorted ( 1 = true)
1, // matrix is "packed" ( 1 = true)
-1, // stype of matrix ( -1 = consider the lower part only )
CHOLMOD_REAL, // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
&c_ // cholmod_common ptr
), deleter);
// copy the data of BCRS matrix to Cholmod Sparse matrix
int* Ap = static_cast<int*>(M->p);
int* Ai = static_cast<int*>(M->i);
Index* Ap = static_cast<Index*>(M->p);
Index* Ai = static_cast<Index*>(M->i);
double* Ax = static_cast<double*>(M->x);
......@@ -287,7 +431,7 @@ public:
// now accumulate
Ap[0] = 0;
for ( int i=0; i<N; i++ )
for ( size_t i=0; i<N; i++ )
{
Ap[i+1] += Ap[i];
}
......@@ -317,14 +461,18 @@ public:
});
// Remove old factor that may be left over from a previous run
if (L_)
CholmodMethod::free_factor(&L_, &c_);
// Now analyse the pattern and optimal row order
L_ = cholmod_analyze(M.get(), &c_);
L_ = CholmodMethod::analyze(M.get(), &c_);
// Do the factorization (this may take some time)
cholmod_factorize(M.get(), L_, &c_);
CholmodMethod::factorize(M.get(), L_, &c_);
}
virtual SolverCategory::Category category() const
SolverCategory::Category category() const override
{
return SolverCategory::Category::sequential;
}
......@@ -332,7 +480,7 @@ public:
/** \brief return a reference to the CHOLMOD common object for advanced option settings
*
* The CHOLMOD common object stores all parameters and options for the solver to run
* and can be modified in several ways, see CHOLMOD Userguide for further information
* and can be modified in several ways, see CHOLMOD Userguide for further information.
*/
cholmod_common& cholmodCommonObject()
{
......@@ -365,7 +513,7 @@ private:
auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
{
const auto deleter = [c](auto* p) {
cholmod_free_dense(&p, c);
CholmodMethod::free_dense(&p, c);
};
return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
}
......@@ -380,34 +528,31 @@ private:
std::vector<std::size_t> subIndices_;
};
struct CholmodCreator{
template<class F> struct isValidBlock : std::false_type{};
template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
template<int k> struct isValidBlock<FieldVector<float,k>> : std::true_type{};
template<class TL, typename M>
std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
typename Dune::TypeListElement<2, TL>::type>>
operator()(TL /*tl*/, const M& mat, const Dune::ParameterTree& /*config*/,
std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
{
using D = typename Dune::TypeListElement<1, TL>::type;
auto solver = std::make_shared<Dune::Cholmod<D>>();
solver->setMatrix(mat);
return solver;
}
// second version with SFINAE to validate the template parameters of Cholmod
template<typename TL, typename M>
std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
typename Dune::TypeListElement<2, TL>::type>>
operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
{
DUNE_THROW(UnsupportedType, "Unsupported Type in Cholmod");
}
};
DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator());
DUNE_REGISTER_SOLVER("cholmod",
[](auto opTraits, const auto& op, const Dune::ParameterTree& config)
-> std::shared_ptr<typename decltype(opTraits)::solver_type>
{
using OpTraits = decltype(opTraits);
using M = typename OpTraits::matrix_type;
using D = typename OpTraits::domain_type;
// works only for sequential operators
if constexpr (OpTraits::isParallel){
if(opTraits.getCommOrThrow(op).communicator().size() > 1)
DUNE_THROW(Dune::InvalidStateException, "CholMod works only for sequential operators.");
}
if constexpr (OpTraits::isAssembled &&
// check whether the Matrix field_type is double or float
(std::is_same_v<typename FieldTraits<D>::field_type, double> ||
std::is_same_v<typename FieldTraits<D>::field_type, float>)){
const auto& A = opTraits.getAssembledOpOrThrow(op);
const M& mat = A->getmat();
auto solver = std::make_shared<Dune::Cholmod<D>>();
solver->setMatrix(mat);
return solver;
}
DUNE_THROW(UnsupportedType,
"Unsupported Type in Cholmod (only double and float supported)");
});
} /* namespace Dune */
......
......@@ -53,7 +53,7 @@ namespace Dune {
}
/*
Register all creators from the registry in the Parameterizedobjectfactory An
Register all creators from the registry in the Parameterizedobjectfactory. An
object of V is passed in the creator and should be used to determine the
template arguments.
*/
......
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_DILU_HH
#define DUNE_ISTL_DILU_HH
#include <cmath>
#include <complex>
#include <map>
#include <vector>
#include <sstream>
#include <dune/common/fmatrix.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include "istlexception.hh"
/** \file
* \brief The diagonal incomplete LU factorization kernels
*/
namespace Dune
{
/** @addtogroup ISTL_Kernel
@{
*/
namespace DILU
{
/*! compute DILU decomposition of A
The preconditioner matrix M has the property
diag(A) = diag(M) = diag((D + L_A) D^{-1} (D + U_A)) = diag(D + L_A D^{-1} U_A)
such that the diagonal matrix D can be constructed:
D_11 = A_11
D_22 = A22 - L_A_{21} D_{11}^{-1} U_A_{12}
and etc
we store the inverse of D to be used when applying the preconditioner.
For more details, see: R. Barrett et al., "Templates for the Solution of Linear Systems:
Building Blocks for Iterative Methods", 1994. Available from: https://www.netlib.org/templates/templates.pdf
*/
template <class M>
void blockDILUDecomposition(M &A, std::vector<typename M::block_type> &Dinv_)
{
auto endi = A.end();
for (auto row = A.begin(); row != endi; ++row)
{
const auto row_i = row.index();
const auto col = row->find(row_i);
// initialise Dinv[i] = A[i, i]
Dinv_[row_i] = *col;
}
for (auto row = A.begin(); row != endi; ++row)
{
const auto row_i = row.index();
for (auto a_ij = row->begin(); a_ij.index() < row_i; ++a_ij)
{
const auto col_j = a_ij.index();
const auto a_ji = A[col_j].find(row_i);
// if A[i, j] != 0 and A[j, i] != 0
if (a_ji != A[col_j].end())
{
// Dinv[i] -= A[i, j] * d[j] * A[j, i]
Dinv_[row_i] -= (*a_ij) * Dinv_[col_j] * (*a_ji);
}
}
// store the inverse
try
{
Impl::asMatrix(Dinv_[row_i]).invert(); // compute inverse of diagonal block
}
catch (Dune::FMatrixError &e)
{
std::ostringstream sstream;
sstream << THROWSPEC(MatrixBlockError)
<< "DILU failed to invert matrix block D[" << row_i << "]" << e.what();
MatrixBlockError ex;
ex.message(sstream.str());
ex.r = row_i;
throw ex;
}
}
}
/*! DILU backsolve
M = (D + L_A) D^-1 (D + U_A) (a LU decomposition of M)
where L_A and U_A are the strictly lower and upper parts of A.
Working with residual d = b - Ax and update v = x_{n+1} - x_n
solving the product M^-1(d) using upper and lower triangular solve:
v = M^{-1} d = (D + U_A)^{-1} D (D + L_A)^{-1} d
define y = (D + L_A)^{-1} d
lower triangular solve: (D + L_A) y = d
upper triangular solve: (D + U_A) v = D y
*/
template <class M, class X, class Y>
void blockDILUBacksolve(const M &A, const std::vector<typename M::block_type> Dinv_, X &v, const Y &d)
{
using dblock = typename Y::block_type;
using vblock = typename X::block_type;
// lower triangular solve: (D + L_A) y = d
auto endi = A.end();
for (auto row = A.begin(); row != endi; ++row)
{
const auto row_i = row.index();
dblock rhsValue(d[row_i]);
auto &&rhs = Impl::asVector(rhsValue);
for (auto a_ij = (*row).begin(); a_ij.index() < row_i; ++a_ij)
{
// if A[i][j] != 0
// rhs -= A[i][j]* y[j], where v_j stores y_j
const auto col_j = a_ij.index();
Impl::asMatrix(*a_ij).mmv(v[col_j], rhs);
}
// y_i = Dinv_i * rhs
// storing y_i in v_i
auto &&vi = Impl::asVector(v[row_i]);
Impl::asMatrix(Dinv_[row_i]).mv(rhs, vi); // (D + L_A)_ii = D_i
}
// upper triangular solve: (D + U_A) v = Dy
auto rendi = A.beforeBegin();
for (auto row = A.beforeEnd(); row != rendi; --row)
{
const auto row_i = row.index();
// rhs = 0
vblock rhs(0.0);
for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij)
{
// if A[i][j] != 0
// rhs += A[i][j]*v[j]
const auto col_j = a_ij.index();
Impl::asMatrix(*a_ij).umv(v[col_j], rhs);
}
// calculate update v = M^-1*d
// v_i = y_i - Dinv_i*rhs
// before update v_i is y_i
auto &&vi = Impl::asVector(v[row_i]);
Impl::asMatrix(Dinv_[row_i]).mmv(rhs, vi);
}
}
} // end namespace DILU
/** @} end documentation */
} // end namespace
#endif
......@@ -204,17 +204,17 @@ namespace Dune
/**
* \brief Wrapper to use a range of ARPACK++ eigenvalue solvers
*
* A class template for performing some eigenvalue algorithms
* provided by the ARPACK++ library which is based on the implicitly
* restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
* Arnoldi/Lanczos process with the implicitily shifted QR technique.
* The method is designed to compute eigenvalue-eigenvector pairs of
* large scale sparse nonsymmetric/symmetric matrices. This class
* template uses the algorithms to compute the dominant (i.e. largest
* magnitude) and least dominant (i.e. smallest magnitude) eigenvalue
* as well as the spectral condition number of square, symmetric
* matrices and to compute the largest and smallest singular value as
* well as the spectral condition number of nonsymmetric matrices.
* A class template for performing some eigenvalue algorithms
* provided by the ARPACK++ library which is based on the implicitly
* restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
* Arnoldi/Lanczos process with the implicitily shifted QR technique.
* The method is designed to compute eigenvalue-eigenvector pairs of
* large scale sparse nonsymmetric/symmetric matrices. This class
* template uses the algorithms to compute the dominant (i.e. largest
* magnitude) and least dominant (i.e. smallest magnitude) eigenvalue
* as well as the spectral condition number of square, symmetric
* matrices and to compute the largest and smallest singular value as
* well as the spectral condition number of nonsymmetric matrices.
*
* \note For a recent version of the ARPACK++ library working with recent
* compiler versions see "http://reuter.mit.edu/software/arpackpatch/"
......@@ -235,8 +235,8 @@ namespace Dune
* respectively singular values shall be considered;
* is assumed to have blocklevel 2.
* \tparam BlockVector Type of the associated vectors; compatible with the
* rows of a BCRSMatrix object (if #rows >= #ncols) or
* its columns (if #rows < #ncols).
* rows of a BCRSMatrix object (if \#rows \>= \#ncols) or
* its columns (if \#rows \< \#ncols).
*
* \author Sebastian Westerheide.
*/
......@@ -256,14 +256,14 @@ namespace Dune
* update iterations allowed; depending on the
* algorithm, c*nIterationsMax iterations may
* be performed, where c is a natural number.
* \param[in] verbosity_level Verbosity setting;
* >= 1: algorithms print a preamble and
* the final result,
* >= 2: algorithms print information about
* the problem solved using ARPACK++,
* >= 3: the final result output includes
* the approximated eigenvector,
* >= 4: sets the ARPACK(++) verbosity mode.
* \param[in] verbosity_level Verbosity setting;<br>
* \>= 1: algorithms print a preamble and
* the final result,<br>
* \>= 2: algorithms print information about
* the problem solved using ARPACK++,<br>
* \>= 3: the final result output includes
* the approximated eigenvector,<br>
* \>= 4: sets the ARPACK(++) verbosity mode.
*/
ArPackPlusPlus_Algorithms (const BCRSMatrix& m,
const unsigned int nIterationsMax = 100000,
......@@ -603,8 +603,8 @@ namespace Dune
* (0 == machine precision).
* \param[out] sigma The approximated largest singular value.
* \param[out] x The associated approximated right-singular
* vector (if #rows >= #ncols) respectively
* left-singular vector (if #rows < #ncols).
* vector (if \#rows \>= \#ncols) respectively
* left-singular vector (if \#rows \< \#ncols).
*/
inline void computeNonSymMax (const Real& epsilon,
BlockVector& x, Real& sigma) const
......@@ -715,8 +715,8 @@ namespace Dune
* (0 == machine precision).
* \param[out] sigma The approximated smallest singular value.
* \param[out] x The associated approximated right-singular
* vector (if #rows >= #ncols) respectively
* left-singular vector (if #rows < #ncols).
* vector (if \#rows \>= \#ncols) respectively
* left-singular vector (if \#rows \< \#ncols).
*/
inline void computeNonSymMin (const Real& epsilon,
BlockVector& x, Real& sigma) const
......