Not getting same result as The inverse-iteration program example

[NickChlanda]

Hello,

I am currently trying to run the inverse-iteration example and am having toruble getting the posted result. My source code is pasted at the bottom of this discussion.

After compiling:

g++ test.cpp -L/usr/local/lib -lginkgo -lginkgo_omp -lginkgo_cuda -lginkgo_hip -lginkgo_dpcpp -lginkgo_device -lginkgo_reference -o test
then running:

./test
I get the result:

`This is Ginkgo 1.10.0 (develop)
running with core module 1.10.0 (develop)
the reference module is 1.10.0 (develop)
the OpenMP module is 1.10.0 (develop)
the CUDA module is 1.10.0 (develop)
the HIP module is not compiled
the DPCPP module is not compiled'

'{ "system_residual": +1.59471978e-14, "eigenvalue": (-3.03897455e+00,-4.55595319e+00), "residual": (-2.49101122e+00,+3.73446053e+00) },
{ "system_residual": +5.32511725e-15, "eigenvalue": (+1.29801689e-02,+1.76338770e+01), "residual": (+1.09858125e-03,-1.49244950e+00) },
{ "system_residual": +5.50869574e-15, "eigenvalue": (+9.91827902e+00,-1.48922780e+01), "residual": (+5.46299627e-01,+8.20267900e-01) },
{ "system_residual": +4.70222503e-15, "eigenvalue": (-1.86260494e+01,+3.54194853e+00), "residual": (-1.97194952e+00,-3.74987928e-01) },
{ "system_residual": +2.67147604e-15, "eigenvalue": (+1.65252354e+01,+9.74791627e+00), "residual": (+4.54423954e-01,-2.68055889e-01) },
{ "system_residual": +3.86242002e-15, "eigenvalue": (-5.60615426e+00,-1.83387674e+01), "residual": (-4.70700469e-01,+1.53974829e+00) },
{ "system_residual": +2.32936948e-15, "eigenvalue": (-8.14431141e+00,+1.74150948e+01), "residual": (-7.14909647e-01,-1.52870128e+00) },
{ "system_residual": +2.09559067e-15, "eigenvalue": (+1.77232613e+01,-7.44163570e+00), "residual": (+3.64196976e-01,+1.52918878e-01) },
{ "system_residual": +1.42776864e-15, "eigenvalue": (-1.81397272e+01,-6.38056937e+00), "residual": (-1.85978092e+00,+6.54169773e-01) },
{ "system_residual": +2.19489239e-15, "eigenvalue": (+9.16599661e+00,+1.69034307e+01), "residual": (+4.87679370e-01,-8.99351680e-01) },
{ "system_residual": +1.41106914e-15, "eigenvalue": (+4.54933804e+00,-1.86832191e+01), "residual": (+2.92337797e-01,+1.20057271e+00) },
{ "system_residual": +2.75150711e-15, "eigenvalue": (-1.59111883e+01,+1.07978967e+01), "residual": (-1.58190660e+00,-1.07353792e+00) },
{ "system_residual": +1.91381942e-15, "eigenvalue": (+1.90424639e+01,+2.67322265e+00), "residual": (+1.38004698e-01,-1.93734007e-02) },
{ "system_residual": +2.18654036e-15, "eigenvalue": (-1.23231791e+01,-1.47614867e+01), "residual": (-1.16094577e+00,+1.39065459e+00) },
{ "system_residual": +3.31640094e-15, "eigenvalue": (-7.70822863e-01,+1.92137451e+01), "residual": (-5.78153307e-02,-1.44112102e+00) },
{ "system_residual": +2.65209348e-15, "eigenvalue": (+1.34660808e+01,-1.37268688e+01), "residual": (+5.42179639e-01,+5.52679648e-01) },
{ "system_residual": +2.68682089e-15, "eigenvalue": (-1.91954298e+01,+1.13917733e+00), "residual": (-1.99561049e+00,-1.18432056e-01) },
{ "system_residual": +2.74704036e-15, "eigenvalue": (+1.49951181e+01,+1.20378028e+01), "residual": (+5.17490939e-01,-4.15432132e-01) },
{ "system_residual": +1.34086489e-13, "eigenvalue": (-3.03793669e+00,-1.89877119e+01), "residual": (-2.40427422e-01,+1.50271947e+00) },
{ "system_residual": +1.14551139e-15, "eigenvalue": (-1.04907481e+01,+1.61154108e+01), "residual": (-9.59187620e-01,-1.47346047e+00) },
`

Here is the Matrix I ran the code with, directly from the examples directory within ginkgo.

MatrixA.txt

which is different from the values from the example in which it converges to ~19.

Source code:

`
#include
#include
#include
#include
#include
#include
#include

#include <ginkgo/ginkgo.hpp>

int main(int argc, char* argv[])
{
using precision = std::complex;
using real_precision = gko::remove_complex;
using vec = gko::matrix::Dense;
using real_vec = gko::matrix::Dense<real_precision>;
using mtx = gko::matrix::Csr;
using solver_type = gko::solver::Bicgstab;

using std::abs;
using std::sqrt;

std::cout << [gko::version_info::get](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1version__info.html#a6daeb8a087cfb57fa055526fc133d8eb)() << std::endl;

std::cout << std::scientific << std::setprecision(8) << std::showpos;

if (argc == 2 && (std::string(argv[1]) == "--help")) {
    std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
    std::exit(-1);
}

const auto executor_string = argc >= 2 ? argv[1] : "reference";
std::map<std::string, std::function<std::shared_ptr<gko::Executor>()>>
    exec_map{
        {"omp", [] { return [gko::OmpExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1OmpExecutor.html#a8930142a1913bcbcac4fb6c917f930d9)(); }},
        {"cuda",
         [] {
             return [gko::CudaExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1CudaExecutor.html#a6c423656c7c0eb53810cb77bb4a66d15)(0,
                                              [gko::OmpExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1OmpExecutor.html#a8930142a1913bcbcac4fb6c917f930d9)());
         }},
        {"hip",
         [] {
             return [gko::HipExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1HipExecutor.html#aa0701fe1303e03a8331da7e95008100d)(0, [gko::OmpExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1OmpExecutor.html#a8930142a1913bcbcac4fb6c917f930d9)());
         }},
        {"dpcpp",
         [] {
             return [gko::DpcppExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1DpcppExecutor.html#ae9b875dfe3af7f3f3f24071bfb2f895e)(0,
                                               [gko::OmpExecutor::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1OmpExecutor.html#a8930142a1913bcbcac4fb6c917f930d9)());
         }},
        {"reference", [] { return gko::ReferenceExecutor::create(); }}};

const auto exec = exec_map.at(executor_string)();  // throws if not valid

auto this_exec = exec->get_master();

auto system_max_iterations = 100u;
auto system_residual_goal = real_precision{1e-16};

auto max_iterations = 20u;
auto residual_goal = real_precision{1e-8};
auto z = precision{20.0, 2.0};

auto A = [share](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a3ce296f73db0ff398bdea6009a3a5c58)(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));

auto [one](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a0059e27f8f4bc348ff65c1e60caf47c8) = [share](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a3ce296f73db0ff398bdea6009a3a5c58)(gko::initialize<vec>({precision{1.0}}, exec));
auto neg_one = [share](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a3ce296f73db0ff398bdea6009a3a5c58)(gko::initialize<vec>({-precision{1.0}}, exec));
auto neg_z = gko::initialize<vec>({-z}, exec);

auto system_matrix = [share](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a3ce296f73db0ff398bdea6009a3a5c58)([gko::Combination<precision>::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1Combination.html)(
    one, A, gko::initialize<vec>({-z}, exec),
    [gko::matrix::Identity<precision>::create](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1matrix_1_1Identity.html)(exec, A->get_size()[0])));

auto [solver](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko_1_1log.html#a15f71654c6c4f4defc95f9359343fa27a7cbf673200498a0d1324f4330825202a) =
    solver_type::build()
        .with_criteria(gko::stop::Iteration::build().with_max_iters(
                           system_max_iterations),
                       [gko::stop::ResidualNorm<precision>::build](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/classgko_1_1stop_1_1ResidualNorm.html)()
                           .with_reduction_factor(system_residual_goal))
        .on(exec)
        ->generate(system_matrix);


auto x = [&] {
    auto work = vec::create(this_exec, [gko::dim<2>](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/structgko_1_1dim.html){A->get_size()[0], 1});
    const auto n = work->get_size()[0];
    for (int i = 0; i < n; ++i) {
        work->get_values()[i] = precision{1.0} / sqrt(n);
    }
    return [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(exec, work);
}();
auto y = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(x);
auto tmp = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(x);
auto norm = gko::initialize<real_vec>({1.0}, exec);
auto inv_norm = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(this_exec, one);
auto g = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(one);

for (auto i = 0u; i < max_iterations; ++i) {
    std::cout << "{ ";
    [solver](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko_1_1log.html#a15f71654c6c4f4defc95f9359343fa27a7cbf673200498a0d1324f4330825202a)->apply(x, y);
    system_matrix->apply(one, y, neg_one, x);
    x->compute_norm2(norm);
    std::cout << "\"system_residual\": "
              << [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(this_exec, norm)->get_values()[0] << ", ";
    x->copy_from(y);
    x->compute_norm2(norm);
    inv_norm->get_values()[0] =
        real_precision{1.0} / [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(this_exec, norm)->get_values()[0];
    x->scale([clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(exec, inv_norm));
    A->apply(x, tmp);
    x->compute_dot(tmp, g);
    auto g_val = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(this_exec, g)->get_values()[0];
    std::cout << "\"eigenvalue\": " << g_val << ", ";
    auto v = gko::initialize<vec>({-g_val}, exec);
    tmp->add_scaled(v, x);
    tmp->compute_norm2(norm);
    auto res_val = [clone](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a1beb80750459e4201aa9d882d2d074c3)(exec->get_master(), norm)->get_values()[0];
    std::cout << "\"residual\": " << res_val / g_val << " }," << std::endl;
    if ([abs](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a99791881a0d56a1c5392ce52b4dce1ac)(res_val) < residual_goal * [abs](https://ginkgo-project.github.io/ginkgo-generated-documentation/doc/develop/namespacegko.html#a99791881a0d56a1c5392ce52b4dce1ac)(g_val)) {
        break;
    }
}

}`

[NickChlanda]

Hello again,

While trying to find an solution I realised I was on Gingko 1.10.0, so, I downloaded the latest version of Ginkgo (1.9.0).

after

cmake -G "Unix Makefiles" .. && cmake --build .
and

cmake --install .

I ran

./inverse-iteration

While in the exampels/inverse-iteration directory and got the same result as printed before.

This is Ginkgo 1.9.0 (main) running with core module 1.9.0 (main) the reference module is 1.9.0 (main) the OpenMP module is 1.9.0 (main) the CUDA module is 1.9.0 (main) the HIP module is not compiled the DPCPP module is not compiled { "system_residual": +1.59471978e-14, "eigenvalue": (-3.03897455e+00,-4.55595319e+00), "residual": (-2.49101122e+00,+3.73446053e+00) }, { "system_residual": +5.32511725e-15, "eigenvalue": (+1.29801689e-02,+1.76338770e+01), "residual": (+1.09858125e-03,-1.49244950e+00) }, { "system_residual": +5.50869574e-15, "eigenvalue": (+9.91827902e+00,-1.48922780e+01), "residual": (+5.46299627e-01,+8.20267900e-01) }, { "system_residual": +4.70222503e-15, "eigenvalue": (-1.86260494e+01,+3.54194853e+00), "residual": (-1.97194952e+00,-3.74987928e-01) }, { "system_residual": +2.67147604e-15, "eigenvalue": (+1.65252354e+01,+9.74791627e+00), "residual": (+4.54423954e-01,-2.68055889e-01) }, { "system_residual": +3.86242002e-15, "eigenvalue": (-5.60615426e+00,-1.83387674e+01), "residual": (-4.70700469e-01,+1.53974829e+00) }, { "system_residual": +2.32936948e-15, "eigenvalue": (-8.14431141e+00,+1.74150948e+01), "residual": (-7.14909647e-01,-1.52870128e+00) }, { "system_residual": +2.09559067e-15, "eigenvalue": (+1.77232613e+01,-7.44163570e+00), "residual": (+3.64196976e-01,+1.52918878e-01) }, { "system_residual": +1.42776864e-15, "eigenvalue": (-1.81397272e+01,-6.38056937e+00), "residual": (-1.85978092e+00,+6.54169773e-01) }, { "system_residual": +2.19489239e-15, "eigenvalue": (+9.16599661e+00,+1.69034307e+01), "residual": (+4.87679370e-01,-8.99351680e-01) }, { "system_residual": +1.41106914e-15, "eigenvalue": (+4.54933804e+00,-1.86832191e+01), "residual": (+2.92337797e-01,+1.20057271e+00) }, { "system_residual": +2.75150711e-15, "eigenvalue": (-1.59111883e+01,+1.07978967e+01), "residual": (-1.58190660e+00,-1.07353792e+00) }, { "system_residual": +1.91381942e-15, "eigenvalue": (+1.90424639e+01,+2.67322265e+00), "residual": (+1.38004698e-01,-1.93734007e-02) }, { "system_residual": +2.18654036e-15, "eigenvalue": (-1.23231791e+01,-1.47614867e+01), "residual": (-1.16094577e+00,+1.39065459e+00) }, { "system_residual": +3.31640094e-15, "eigenvalue": (-7.70822863e-01,+1.92137451e+01), "residual": (-5.78153307e-02,-1.44112102e+00) }, { "system_residual": +2.65209348e-15, "eigenvalue": (+1.34660808e+01,-1.37268688e+01), "residual": (+5.42179639e-01,+5.52679648e-01) }, { "system_residual": +2.68682089e-15, "eigenvalue": (-1.91954298e+01,+1.13917733e+00), "residual": (-1.99561049e+00,-1.18432056e-01) }, { "system_residual": +2.74704036e-15, "eigenvalue": (+1.49951181e+01,+1.20378028e+01), "residual": (+5.17490939e-01,-4.15432132e-01) }, { "system_residual": +1.34086489e-13, "eigenvalue": (-3.03793669e+00,-1.89877119e+01), "residual": (-2.40427422e-01,+1.50271947e+00) }, { "system_residual": +1.14551139e-15, "eigenvalue": (-1.04907481e+01,+1.61154108e+01), "residual": (-9.59187620e-01,-1.47346047e+00) },

[NickChlanda]

Hello,

I have found the problem in the code.

from the example code I have replaced

x->compute_dot(tmp, g); auto g_val = clone(this_exec, g) ->get_values() [0];

with

precision g_val{0.0}; for (gko::size_type j = 0; j < x->get_size()[0]; ++j) { g_val += std::conj(x->get_values()[j]) * tmp->get_values()[j]; }

Which then gives the expected output

This is Ginkgo 1.10.0 (develop) running with core module 1.10.0 (develop) the reference module is 1.10.0 (develop) the OpenMP module is 1.10.0 (develop) the CUDA module is 1.10.0 (develop) the HIP module is not compiled the DPCPP module is not compiled

{ "system_residual": +1.59471978e-14, "eigenvalue": (+2.03741410e+01,+3.52582546e-16), "residual": (+2.92231055e-01,-5.05717369e-18) }, { "system_residual": +5.32511725e-15, "eigenvalue": (+1.94878474e+01,+1.72388145e-16), "residual": (+7.94370276e-02,-7.02694432e-19) }, { "system_residual": +5.50869574e-15, "eigenvalue": (+1.93282121e+01,-1.35974578e-15), "residual": (+4.11149623e-02,+2.89245048e-18) }, { "system_residual": +4.70222503e-15, "eigenvalue": (+1.92638912e+01,+3.30777271e-16), "residual": (+2.34717040e-02,-4.03028969e-19) }, { "system_residual": +2.67147604e-15, "eigenvalue": (+1.92409166e+01,+2.46648900e-16), "residual": (+1.34709547e-02,-1.72683882e-19) }, { "system_residual": +3.86242002e-15, "eigenvalue": (+1.92331106e+01,-1.34483739e-15), "residual": (+7.72060707e-03,+5.39848246e-19) }, { "system_residual": +2.32936948e-15, "eigenvalue": (+1.92305014e+01,-9.81687892e-17), "residual": (+4.42106625e-03,+2.25688718e-20) }, { "system_residual": +2.09559067e-15, "eigenvalue": (+1.92296339e+01,+6.81443618e-16), "residual": (+2.53081312e-03,-8.96848298e-20) }, { "system_residual": +1.42776864e-15, "eigenvalue": (+1.92293461e+01,+8.76538281e-17), "residual": (+1.44862114e-03,-6.60330246e-21) }, { "system_residual": +2.19489239e-15, "eigenvalue": (+1.92292506e+01,-1.40989206e-15), "residual": (+8.29183451e-04,+6.07958773e-20) }, { "system_residual": +1.41106914e-15, "eigenvalue": (+1.92292190e+01,+4.60500050e-17), "residual": (+4.74636702e-04,-1.13665680e-21) }, { "system_residual": +2.75150711e-15, "eigenvalue": (+1.92292085e+01,+1.48269835e-16), "residual": (+2.71701077e-04,-2.09499387e-21) }, { "system_residual": +1.91381942e-15, "eigenvalue": (+1.92292051e+01,+1.90322983e-16), "residual": (+1.55539937e-04,-1.53947210e-21) }, { "system_residual": +2.18654036e-15, "eigenvalue": (+1.92292039e+01,-2.96720512e-16), "residual": (+8.90457139e-05,+1.37403971e-21) }, { "system_residual": +3.31640094e-15, "eigenvalue": (+1.92292035e+01,+1.55634681e-16), "residual": (+5.09805252e-05,-4.12619158e-22) }, { "system_residual": +2.65209348e-15, "eigenvalue": (+1.92292034e+01,+4.22300558e-16), "residual": (+2.91887365e-05,-6.41026019e-22) }, { "system_residual": +2.68682089e-15, "eigenvalue": (+1.92292034e+01,+3.62717110e-17), "residual": (+1.67126561e-05,-3.15247918e-23) }, { "system_residual": +2.74704036e-15, "eigenvalue": (+1.92292034e+01,-2.09656854e-16), "residual": (+9.56961199e-06,+1.04337902e-22) }, { "system_residual": +1.34086489e-13, "eigenvalue": (+1.92292033e+01,-7.82886719e-16), "residual": (+5.47975753e-06,+2.23099695e-22) }, { "system_residual": +1.14551139e-15, "eigenvalue": (+1.92292033e+01,-4.98143888e-16), "residual": (+3.13794996e-06,+8.12904499e-23) },

So, something might be wrong with compute_dot()? or how it was being used? I am not entirely sure.

[nbeams]

Ah, looks like the example should call compute_conj_dot rather than compute_dot (which is, of course, what you are doing with your code, but compute_conj_dot will work for all backends). Thanks for catching this!

[nbeams]

Hi @NickChlanda, just wanted to say that the fix for this has just been merged (#1878). Thanks again for bringing it to our attention!