Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mixed Precision ISAI #719

Closed
wants to merge 5 commits into from
Closed

Mixed Precision ISAI #719

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c697171
backup
fritzgoebel Mar 9, 2021
3960d29
Make ISAI support reduced storage precision
fritzgoebel Mar 10, 2021
5f33601
get rid of extra set of types
fritzgoebel Mar 11, 2021
414eb8f
Improve format conversion code in ISAI
Mar 23, 2021
334b5ca
Review comments
fritzgoebel Jun 6, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 49 additions & 26 deletions core/preconditioner/isai.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,9 @@ std::shared_ptr<Csr> extend_sparsity(std::shared_ptr<const Executor> &exec,
}


template <isai_type IsaiType, typename ValueType, typename IndexType>
void Isai<IsaiType, ValueType, IndexType>::generate_inverse(
template <isai_type IsaiType, typename ValueType, typename IndexType,
typename StorageType>
void Isai<IsaiType, ValueType, IndexType, StorageType>::generate_inverse(
std::shared_ptr<const LinOp> input, bool skip_sorting, int power,
IndexType excess_limit)
{
Expand Down Expand Up @@ -265,58 +266,80 @@ void Isai<IsaiType, ValueType, IndexType>::generate_inverse(
}


template <isai_type IsaiType, typename ValueType, typename IndexType>
std::unique_ptr<LinOp> Isai<IsaiType, ValueType, IndexType>::transpose() const
template <isai_type IsaiType, typename ValueType, typename IndexType,
typename StorageType>
std::unique_ptr<LinOp>
Isai<IsaiType, ValueType, IndexType, StorageType>::transpose() const
{
auto is_spd = IsaiType == isai_type::spd;
if (is_spd) {
if (IsaiType == isai_type::spd) {
return this->clone();
}

std::unique_ptr<transposed_type> transp{
new transposed_type{this->get_executor()}};
transp->set_size(gko::transpose(this->get_size()));

transp->approximate_inverse_ =
share(as<Csr>(this->get_approximate_inverse())->transpose());
convert_csr_to_ell(this->get_approximate_inverse()->transpose().get());

return std::move(transp);
}


template <isai_type IsaiType, typename ValueType, typename IndexType>
std::unique_ptr<LinOp> Isai<IsaiType, ValueType, IndexType>::conj_transpose()
const
template <isai_type IsaiType, typename ValueType, typename IndexType,
typename StorageType>
std::unique_ptr<LinOp>
Isai<IsaiType, ValueType, IndexType, StorageType>::conj_transpose() const
{
auto is_spd = IsaiType == isai_type::spd;
if (is_spd) {
if (IsaiType == isai_type::spd) {
return this->clone();
}

std::unique_ptr<transposed_type> transp{
new transposed_type{this->get_executor()}};
transp->set_size(gko::transpose(this->get_size()));
transp->approximate_inverse_ =
share(as<Csr>(this->get_approximate_inverse())->conj_transpose());

transp->approximate_inverse_ = convert_csr_to_ell(
this->get_approximate_inverse()->conj_transpose().get());

return std::move(transp);
}


#define GKO_DECLARE_LOWER_ISAI(ValueType, IndexType) \
class Isai<isai_type::lower, ValueType, IndexType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI);
#define GKO_DECLARE_LOWER_ISAI1(ValueType, IndexType) \
class Isai<isai_type::lower, ValueType, IndexType, ValueType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI1);

#define GKO_DECLARE_UPPER_ISAI1(ValueType, IndexType) \
class Isai<isai_type::upper, ValueType, IndexType, ValueType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI1);

#define GKO_DECLARE_GENERAL_ISAI1(ValueType, IndexType) \
class Isai<isai_type::general, ValueType, IndexType, ValueType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI1);

#define GKO_DECLARE_SPD_ISAI1(ValueType, IndexType) \
class Isai<isai_type::spd, ValueType, IndexType, ValueType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI1);

#define GKO_DECLARE_LOWER_ISAI2(ValueType, IndexType) \
class Isai<isai_type::lower, ValueType, IndexType, \
next_precision<ValueType>>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI2);

#define GKO_DECLARE_UPPER_ISAI(ValueType, IndexType) \
class Isai<isai_type::upper, ValueType, IndexType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI);
#define GKO_DECLARE_UPPER_ISAI2(ValueType, IndexType) \
class Isai<isai_type::upper, ValueType, IndexType, \
next_precision<ValueType>>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI2);

#define GKO_DECLARE_GENERAL_ISAI(ValueType, IndexType) \
class Isai<isai_type::general, ValueType, IndexType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI);
#define GKO_DECLARE_GENERAL_ISAI2(ValueType, IndexType) \
class Isai<isai_type::general, ValueType, IndexType, \
next_precision<ValueType>>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI2);

#define GKO_DECLARE_SPD_ISAI(ValueType, IndexType) \
class Isai<isai_type::spd, ValueType, IndexType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI);
#define GKO_DECLARE_SPD_ISAI2(ValueType, IndexType) \
class Isai<isai_type::spd, ValueType, IndexType, next_precision<ValueType>>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI2);


} // namespace preconditioner
Expand Down
138 changes: 111 additions & 27 deletions include/ginkgo/core/preconditioner/isai.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <ginkgo/core/base/lin_op.hpp>
#include <ginkgo/core/matrix/csr.hpp>
#include <ginkgo/core/matrix/dense.hpp>
#include <ginkgo/core/matrix/ell.hpp>


namespace gko {
Expand Down Expand Up @@ -102,30 +103,31 @@ enum struct isai_type { lower, upper, general, spd };
* @ingroup precond
* @ingroup LinOp
*/
template <isai_type IsaiType, typename ValueType, typename IndexType>
class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
public Transposable {
template <isai_type IsaiType, typename ValueType, typename IndexType,
typename StorageType = ValueType>
class Isai
: public EnableLinOp<Isai<IsaiType, ValueType, IndexType, StorageType>>,
public Transposable {
friend class EnableLinOp<Isai>;
friend class EnablePolymorphicObject<Isai, LinOp>;
friend class Isai<isai_type::general, ValueType, IndexType>;
friend class Isai<isai_type::lower, ValueType, IndexType>;
friend class Isai<isai_type::upper, ValueType, IndexType>;
friend class Isai<isai_type::spd, ValueType, IndexType>;
friend class Isai<isai_type::general, ValueType, IndexType, StorageType>;
friend class Isai<isai_type::lower, ValueType, IndexType, StorageType>;
friend class Isai<isai_type::upper, ValueType, IndexType, StorageType>;
friend class Isai<isai_type::spd, ValueType, IndexType, StorageType>;

public:
using value_type = ValueType;
using index_type = IndexType;
using transposed_type =
Isai<IsaiType == isai_type::general
? isai_type::general
: IsaiType == isai_type::spd
? isai_type::spd
: IsaiType == isai_type::lower ? isai_type::upper
: isai_type::lower,
ValueType, IndexType>;
Isai<IsaiType == isai_type::lower
? isai_type::upper
: IsaiType == isai_type::upper ? isai_type::lower : IsaiType,
Comment on lines +123 to +125
Copy link
Member

@upsj upsj Jun 10, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could we move this into a

constexpr isai_type transpose(isai_type type) {
    return type == isai_type::lower ? isai_type::upper : type == isai_type::upper ? isai_type::lower : type;
}

function?

Suggested change
Isai<IsaiType == isai_type::lower
? isai_type::upper
: IsaiType == isai_type::upper ? isai_type::lower : IsaiType,
Isai<transpose(IsaiType),

ValueType, IndexType, StorageType>;

using Comp = Composition<ValueType>;
using Csr = matrix::Csr<ValueType, IndexType>;
using Dense = matrix::Dense<ValueType>;
using Ell = matrix::Ell<StorageType, IndexType>;
static constexpr isai_type type{IsaiType};

/**
Expand All @@ -139,8 +141,19 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
Comp, Csr>::type>
get_approximate_inverse() const
{
auto exec = this->get_executor();
std::shared_ptr<const LinOp> return_inv;
if (IsaiType == isai_type::spd) {
auto ops = as<Comp>(approximate_inverse_)->get_operators();
auto csr_inv = convert_ell_to_csr(ops[1].get());
auto csr_inv_transp = convert_ell_to_csr(ops[0].get());
return_inv =
share(Composition<ValueType>::create(csr_inv_transp, csr_inv));
} else {
return_inv = convert_ell_to_csr(approximate_inverse_.get());
}
return as<typename std::conditional<IsaiType == isai_type::spd, Comp,
Csr>::type>(approximate_inverse_);
Csr>::type>(return_inv);
}

GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory)
Expand Down Expand Up @@ -186,6 +199,15 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
*/
std::shared_ptr<LinOpFactory> GKO_FACTORY_PARAMETER_SCALAR(
excess_solver_factory, nullptr);

/**
* @brief Parameter for reducing preconditioner storage precision.
*
* Reduces the preconditioner storage precision up to reduce_precision
* times. For value types double and complex<double>, storage precision
* can be reduced up to two times, for float and complex<float> once.
*/
int GKO_FACTORY_PARAMETER_SCALAR(reduce_precision, 0);
};

GKO_ENABLE_LIN_OP_FACTORY(Isai, parameters, Factory);
Expand All @@ -211,15 +233,21 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
: EnableLinOp<Isai>(factory->get_executor(), system_matrix->get_size()),
parameters_{factory->get_parameters()}
{
const auto exec = this->get_executor();
const auto skip_sorting = parameters_.skip_sorting;
const auto power = parameters_.sparsity_power;
const auto excess_limit = parameters_.excess_limit;
generate_inverse(system_matrix, skip_sorting, power, excess_limit);

auto csr_inv = as<Csr>(approximate_inverse_);
auto ell_inv = convert_csr_to_ell(approximate_inverse_.get());
if (IsaiType == isai_type::spd) {
auto inv = share(as<Csr>(approximate_inverse_));
auto inv_transp = share(inv->conj_transpose());
auto csr_inv_transp = csr_inv->conj_transpose();
auto ell_inv_transp = convert_csr_to_ell(csr_inv_transp.get());
approximate_inverse_ =
Composition<ValueType>::create(inv_transp, inv);
Composition<ValueType>::create(ell_inv_transp, ell_inv);
} else {
approximate_inverse_ = ell_inv;
}
}

Expand All @@ -235,6 +263,59 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
}

private:
/**
* Converts a given matrix in ELL format to a matrix in CSR format.
* @note ELL format uses StorageType as the storage
*
* @param in Matrix in ELL format
*
* @returns the given matrix converted to the CSR format.
*/
static std::shared_ptr<Csr> convert_ell_to_csr(const LinOp *in)
{
using EllValue = matrix::Ell<value_type, index_type>;
auto ell = as<const Ell>(in);
auto out = Csr::create(ell->get_executor());
if (std::is_same<value_type, StorageType>::value) {
// This cast is necessary for the compiler
as<EllValue>(ell)->convert_to(out.get());
}
// In-between step required to convert StorageType to value_type
else {
auto tmp = EllValue::create(ell->get_executor());
ell->convert_to(tmp.get());
tmp->convert_to(out.get());
}
return std::move(out);
}

/**
* Converts a given matrix in CSR format to a matrix in ELL format
* @note ELL format uses StorageType as the storage
*
* @param in Matrix in CSR format
*
* @returns the given matrix converted to the ELL format.
*/
static std::shared_ptr<Ell> convert_csr_to_ell(const LinOp *in)
{
using CsrStorage = matrix::Csr<StorageType, index_type>;

auto csr = as<const Csr>(in);
auto out = Ell::create(in->get_executor());
if (std::is_same<value_type, StorageType>::value) {
// This cast is necessary for the compiler
as<CsrStorage>(csr)->convert_to(out.get());
}
// In-between step required to convert value_type to StorageType
else {
auto tmp = CsrStorage::create(csr->get_executor());
csr->convert_to(tmp.get());
tmp->convert_to(out.get());
}
return std::move(out);
}

/**
* Generates the approximate inverse for a triangular matrix and
* stores the result in `approximate_inverse_`.
Expand All @@ -249,22 +330,25 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
bool skip_sorting, int power,
index_type excess_limit);

private:
std::shared_ptr<LinOp> approximate_inverse_;
};


template <typename ValueType = default_precision, typename IndexType = int32>
using LowerIsai = Isai<isai_type::lower, ValueType, IndexType>;
template <typename ValueType = default_precision, typename IndexType = int32,
typename StorageType = ValueType>
using LowerIsai = Isai<isai_type::lower, ValueType, IndexType, StorageType>;

template <typename ValueType = default_precision, typename IndexType = int32>
using UpperIsai = Isai<isai_type::upper, ValueType, IndexType>;
template <typename ValueType = default_precision, typename IndexType = int32,
typename StorageType = ValueType>
using UpperIsai = Isai<isai_type::upper, ValueType, IndexType, StorageType>;

template <typename ValueType = default_precision, typename IndexType = int32>
using GeneralIsai = Isai<isai_type::general, ValueType, IndexType>;
template <typename ValueType = default_precision, typename IndexType = int32,
typename StorageType = ValueType>
using GeneralIsai = Isai<isai_type::general, ValueType, IndexType, StorageType>;

template <typename ValueType = default_precision, typename IndexType = int32>
using SpdIsai = Isai<isai_type::spd, ValueType, IndexType>;
template <typename ValueType = default_precision, typename IndexType = int32,
typename StorageType = ValueType>
using SpdIsai = Isai<isai_type::spd, ValueType, IndexType, StorageType>;


} // namespace preconditioner
Expand Down
Loading