[petsc-users] Multiplication of partitioned with non-partitioned (sparse) PETSc matrices

[Pierre Jolivet]



> On 23 Aug 2023, at 5:35 PM, Thanasis Boutsikakis <thanasis.boutsikakis at corintis.com> wrote:
> Hi all,
> 
> I am trying to multiply two Petsc matrices as C = A * B, where A is a tall matrix and B is a relatively small matrix.
> 
> I have taken the decision to create A as (row-)partitioned matrix and B as a non-partitioned matrix that it is entirely shared by all procs (to avoid unnecessary communication).
> 
> Here is my code:
> 
> import numpy as np
> from firedrake import COMM_WORLD
> from firedrake.petsc import PETSc
> from numpy.testing import assert_array_almost_equal
> 
> nproc = COMM_WORLD.size
> rank = COMM_WORLD.rank
> 
> def create_petsc_matrix_non_partitioned(input_array):
>     """Building a mpi non-partitioned petsc matrix from an array
> 
>     Args:
>         input_array (np array): Input array
>         sparse (bool, optional): Toggle for sparese or dense. Defaults to True.
> 
>     Returns:
>         mpi mat: PETSc matrix
>     """
>     assert len(input_array.shape) == 2
> 
>     m, n = input_array.shape
> 
>     matrix = PETSc.Mat().createAIJ(size=((m, n), (m, n)), comm=COMM_WORLD)
> 
>     # Set the values of the matrix
>     matrix.setValues(range(m), range(n), input_array[:, :], addv=False)
> 
>     # Assembly the matrix to compute the final structure
>     matrix.assemblyBegin()
>     matrix.assemblyEnd()
> 
>     return matrix
> 
> 
> def create_petsc_matrix(input_array, partition_like=None):
>     """Create a PETSc matrix from an input_array
> 
>     Args:
>         input_array (np array): Input array
>         partition_like (petsc mat, optional): Petsc matrix. Defaults to None.
>         sparse (bool, optional): Toggle for sparese or dense. Defaults to True.
> 
>     Returns:
>         petsc mat: PETSc matrix
>     """
>     # Check if input_array is 1D and reshape if necessary
>     assert len(input_array.shape) == 2, "Input array should be 2-dimensional"
>     global_rows, global_cols = input_array.shape
> 
>     comm = COMM_WORLD
>     if partition_like is not None:
>         local_rows_start, local_rows_end = partition_like.getOwnershipRange()
>         local_rows = local_rows_end - local_rows_start
> 
>         # No parallelization in the columns, set local_cols = None to parallelize
>         size = ((local_rows, global_rows), (global_cols, global_cols))
>     else:
>         size = ((None, global_rows), (global_cols, global_cols))
> 
>     matrix = PETSc.Mat().createAIJ(size=size, comm=comm)
>     matrix.setUp()
> 
>     local_rows_start, local_rows_end = matrix.getOwnershipRange()
> 
>     for counter, i in enumerate(range(local_rows_start, local_rows_end)):
>         # Calculate the correct row in the array for the current process
>         row_in_array = counter + local_rows_start
>         matrix.setValues(
>             i, range(global_cols), input_array[row_in_array, :], addv=False
>         )
> 
>     # Assembly the matrix to compute the final structure
>     matrix.assemblyBegin()
>     matrix.assemblyEnd()
> 
>     return matrix
> 
> 
> m, k = 10, 3
> # Generate the random numpy matrices
> np.random.seed(0)  # sets the seed to 0
> A_np = np.random.randint(low=0, high=6, size=(m, k))
> B_np = np.random.randint(low=0, high=6, size=(k, k))
> 
> 
> A = create_petsc_matrix(A_np)
> 
> B = create_petsc_matrix_non_partitioned(B_np)
> 
> # Now perform the multiplication
> C = A * B
> 
> The problem with this is that there is a mismatch between the local rows of A (depend on the partitioning) and the global rows of B (3 for all procs), so that the multiplication cannot happen in parallel. Here is the error:
> 
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
> [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
> [0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
> [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
> [0]PETSC ERROR: to get more information on the crash.
> [0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
> [unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59
> :
> system msg for write_line failure : Bad file descriptor
> 
> 
> Is there a standard way to achieve this?

Your B is duplicated by all processes?
If so, then, call https://petsc.org/main/manualpages/Mat/MatMPIAIJGetLocalMat/, do a sequential product with B on COMM_SELF, not COMM_WORLD, and use https://petsc.org/main/manualpages/Mat/MatCreateMPIMatConcatenateSeqMat/ with the output.

Thanks,
Pierre

> Thanks,
> Thanos
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230823/c0a607c6/attachment-0001.html>