Skip to content

Commit

Permalink
review updates:
Browse files Browse the repository at this point in the history 
- documentation
- formatting
- fix reduction type issue

Co-authored-by: Yu-Hsiang M. Tsai <yhmtsai@gmail.com>
  • Loading branch information
@MarcelKoch @yhmtsai
MarcelKoch and yhmtsai committed Apr 24, 2023
1 parent e3e2d73 commit 6a6d4b7
Show file tree
Hide file tree
Showing 4 changed files with 26 additions and 25 deletions.
40 changes: 21 additions & 19 deletions core/test/utils/matrix_generator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -509,25 +509,25 @@ std::unique_ptr<MatrixType> generate_random_band_matrix(
* diag, upper]
* @param exec the executor for the resulting matrix
*
* @ return a tridiagonal Toeplitz matrix.
* @return a tridiagonal Toeplitz matrix.
*/
template <typename ValueType, typename IndexType>
gko::matrix_data<ValueType, IndexType> generate_tridiag_matrix_data(
gko::size_type size, std::array<ValueType, 3> coeffs,
std::shared_ptr<const gko::Executor> exec)
{
auto a = coeffs[1];
auto b = coeffs[2];
auto c = coeffs[0];
auto lower = coeffs[0];
auto diag = coeffs[1];
auto upper = coeffs[2];

gko::matrix_data<ValueType, IndexType> md{gko::dim<2>{size, size}};
for (int i = 0; i < size; ++i) {
for (size_type i = 0; i < size; ++i) {
if (i > 0) {
md.nonzeros.emplace_back(i, i - 1, c);
md.nonzeros.emplace_back(i, i - 1, lower);
}
md.nonzeros.emplace_back(i, i, a);
md.nonzeros.emplace_back(i, i, diag);
if (i < size - 1) {
md.nonzeros.emplace_back(i, i + 1, b);
md.nonzeros.emplace_back(i, i + 1, upper);
}
}
return md;
Expand Down Expand Up @@ -561,41 +561,43 @@ std::unique_ptr<MatrixType> generate_tridiag_matrix(
* diag, upper]
* @param exec the executor for the resulting matrix
*
* @ return a matrix (possible dense) that is the inverse of the matrix
* @return a matrix (possible dense) that is the inverse of the matrix
* generated from generate_tridiag_matrix_data with the same inputs
*/
template <typename ValueType, typename IndexType>
gko::matrix_data<ValueType, IndexType> generate_tridiag_inverse_matrix_data(
gko::size_type size, std::array<ValueType, 3> coeffs,
std::shared_ptr<const gko::Executor> exec)
{
auto a = coeffs[1];
auto b = coeffs[2];
auto c = coeffs[0];
auto lower = coeffs[0];
auto diag = coeffs[1];
auto upper = coeffs[2];

std::vector<ValueType> alpha(size + 1);
auto beta = std::make_reverse_iterator(alpha.end());

alpha[0] = 1;
alpha[1] = a;
for (int i = 2; i < alpha.size(); ++i) {
alpha[i] = a * alpha[i - 1] - b * c * alpha[i - 2];
alpha[1] = diag;
for (std::size_t i = 2; i < alpha.size(); ++i) {
alpha[i] = diag * alpha[i - 1] - upper * lower * alpha[i - 2];
}

gko::matrix_data<ValueType, IndexType> md{gko::dim<2>{size, size}};
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
for (size_type i = 0; i < size; ++i) {
for (size_type j = 0; j < size; ++j) {
if (i == j) {
md.nonzeros.emplace_back(i, j,
alpha[i] * beta[j + 1] / alpha.back());
} else if (i < j) {
auto sign = static_cast<ValueType>((i + j) % 2 ? -1 : 1);
auto val = sign * static_cast<ValueType>(std::pow(b, j - i)) *
auto val = sign *
static_cast<ValueType>(std::pow(upper, j - i)) *
alpha[i] * beta[j + 1] / alpha.back();
md.nonzeros.emplace_back(i, j, val);
} else {
auto sign = static_cast<ValueType>((i + j) % 2 ? -1 : 1);
auto val = sign * static_cast<ValueType>(std::pow(c, i - j)) *
auto val = sign *
static_cast<ValueType>(std::pow(lower, i - j)) *
alpha[j] * beta[i + 1] / alpha.back();
md.nonzeros.emplace_back(i, j, val);
}
Expand Down
7 changes: 3 additions & 4 deletions include/ginkgo/core/preconditioner/isai.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,9 +213,7 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
std::shared_ptr<LinOpFactory> GKO_FACTORY_PARAMETER_SCALAR(
excess_solver_factory, nullptr);

gko::remove_complex<value_type> GKO_FACTORY_PARAMETER_SCALAR(
excess_solver_reduction,
static_cast<gko::remove_complex<value_type>>(1e-6));
double GKO_FACTORY_PARAMETER_SCALAR(excess_solver_reduction, 1e-6);
};

GKO_ENABLE_LIN_OP_FACTORY(Isai, parameters, Factory);
Expand Down Expand Up @@ -245,7 +243,8 @@ class Isai : public EnableLinOp<Isai<IsaiType, ValueType, IndexType>>,
const auto power = parameters_.sparsity_power;
const auto excess_limit = parameters_.excess_limit;
generate_inverse(system_matrix, skip_sorting, power, excess_limit,
parameters_.excess_solver_reduction);
static_cast<remove_complex<value_type>>(
parameters_.excess_solver_reduction));
if (IsaiType == isai_type::spd) {
auto inv = share(as<Csr>(approximate_inverse_));
auto inv_transp = share(inv->conj_transpose());
Expand Down
1 change: 1 addition & 0 deletions reference/test/preconditioner/isai_kernels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1622,4 +1622,5 @@ TYPED_TEST(Isai, IsExactInverseOnFullSparsitySetLarge)
5 * r<value_type>::value);
}


} // namespace
3 changes: 1 addition & 2 deletions test/preconditioner/isai_kernels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -663,7 +663,6 @@ TEST_F(Isai, IsaiGenerateGeneralWithSparsityPower8IsEquivalentToRef)

auto isai =
Isai::build().with_sparsity_power(8).on(ref)->generate(mtx->clone());

auto d_isai =
Isai::build().with_sparsity_power(8).on(exec)->generate(d_mtx->clone());

Expand All @@ -672,6 +671,7 @@ TEST_F(Isai, IsaiGenerateGeneralWithSparsityPower8IsEquivalentToRef)
r<value_type>::value);
}


TEST_F(Isai, IsaiGenerateGeneralSparsityPowerNIsEquivalentToRef)
{
using Isai =
Expand All @@ -684,7 +684,6 @@ TEST_F(Isai, IsaiGenerateGeneralSparsityPowerNIsEquivalentToRef)
.with_excess_solver_reduction(r<value_type>::value)
.on(ref)
->generate(mtx->clone());

auto d_isai =
Isai::build()
.with_sparsity_power(static_cast<int>(d_mtx->get_size()[0]))
Expand Down

0 comments on commit 6a6d4b7

Please sign in to comment.