STRUMPACK icon indicating copy to clipboard operation
STRUMPACK copied to clipboard

construct_partially_matrix_free for MPI

Open AlexGKim opened this issue 1 year ago • 2 comments

I would like to check whether MPI StructuredMatrix construction using construct_partially_matrix_free is implemented.

I added the following test to examples/dense/testStructuredMPI.cc , which compiles but terminates with an error.

` if (world.is_root()) cout << endl << endl << "====================================" << endl << " Compression, partially matrix-free " << endl << "====================================" << endl; for (auto type : types) { options.set_type(type); try { auto Toeplitz_block = [&Toeplitz](const std::vectorstd::size_t& I, const std::vectorstd::size_t& J, DistributedMatrix& B) { for (std::size_t j=0; j<J.size(); j++) for (std::size_t i=0; i<I.size(); i++) B(i, j) = Toeplitz(I[i], J[j]); };

        // Construct a StructuredMatrix using both a (fast)
  // matrix-vector multiplication and an element extraction
  // routine. This is mainly usefull for HSS construction, which
  // requires element extraction for the diagonal blocks and
  // random projection with the matrix-vector multiplication for
  // the off-diagonal compression.
  auto H = structured::construct_partially_matrix_free<double>
    (&grid, n, n, Tmult2d, Toeplitz_block, options);

   // print_info(world, H.get(), options);
   // check_accuracy(A2d, H.get());
   // factor_and_solve(nrhs, A2d, H.get());
   // preconditioned_solve(world, Tmult1d, H.get());
   // test_shift(nrhs, A2d, H.get());
} catch (std::exception& e) {
    if (world.is_root())
      cout << get_name(type) << " failed: " << e.what() << endl;
  }
}`

AlexGKim avatar Jul 03 '23 16:07 AlexGKim

Try replacing

B(i, j) = Toeplitz(I[i], J[j]);

with

B.global(i, j, Toeplitz(I[i], J[j]));

the global routine treats i and j as global indices in the distributed matrix B, and it will only assign the value if i, j is local to B.

Or equivalently, use:

if (B.is_local(i, j))
  B.global(i, j) = Toeplitz(I[i], J[j]);

then the check is done separately, so you don't need to compute the value if it's not locally required.

pghysels avatar Jul 03 '23 18:07 pghysels

This works, thanks!

AlexGKim avatar Jul 03 '23 21:07 AlexGKim