Hi all, Revisiting my code and the proposed solution from Pierre, I realized this works only in sequential. The reason is that PETSc partitions those matrices only row-wise, which leads to an error due to the mismatch between number of columns of A (non-partitioned) and the number of rows of Phi (partitioned).
"""Experimenting with PETSc mat-mat multiplication""" import time import numpy as np from colorama import Fore from firedrake import COMM_SELF, COMM_WORLD from firedrake.petsc import PETSc from mpi4py import MPI from numpy.testing import assert_array_almost_equal from utilities import Print nproc = COMM_WORLD.size rank = COMM_WORLD.rank def create_petsc_matrix(input_array, sparse=True): """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 mpi 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 size = ((None, global_rows), (global_cols, global_cols)) # Create a sparse or dense matrix based on the 'sparse' argument if sparse: matrix = PETSc.Mat().createAIJ(size=size, comm=COMM_WORLD) else: matrix = PETSc.Mat().createDense(size=size, comm=COMM_WORLD) 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 # -------------------------------------------- # EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc matrix Phi # A' = Phi.T * A * Phi # [k x k] <- [k x m] x [m x m] x [m x k] # -------------------------------------------- m, k = 100, 7 # 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, m)) Phi_np = np.random.randint(low=0, high=6, size=(m, k)) # -------------------------------------------- # TEST: Galerking projection of numpy matrices A_np and Phi_np # -------------------------------------------- Aprime_np = Phi_np.T @ A_np @ Phi_np Print(f"MATRIX Aprime_np [{Aprime_np.shape[0]}x{Aprime_np.shape[1]}]") Print(f"{Aprime_np}") # Create A as an mpi matrix distributed on each process A = create_petsc_matrix(A_np, sparse=False) # Create Phi as an mpi matrix distributed on each process Phi = create_petsc_matrix(Phi_np, sparse=False) # Create an empty PETSc matrix object to store the result of the PtAP operation. # This will hold the result A' = Phi.T * A * Phi after the computation. A_prime = create_petsc_matrix(np.zeros((k, k)), sparse=False) # Perform the PtAP (Phi Transpose times A times Phi) operation. # In mathematical terms, this operation is A' = Phi.T * A * Phi. # A_prime will store the result of the operation. A_prime = A.ptap(Phi) Here is the error MATRIX mpiaij A [100x100] Assembled Partitioning for A: Rank 0: Rows [0, 34) Rank 1: Rows [34, 67) Rank 2: Rows [67, 100) MATRIX mpiaij Phi [100x7] Assembled Partitioning for Phi: Rank 0: Rows [0, 34) Rank 1: Rows [34, 67) Rank 2: Rows [67, 100) Traceback (most recent call last): File "/Users/boutsitron/work/galerkin_projection.py", line 87, in <module> A_prime = A.ptap(Phi) ^^^^^^^^^^^ File "petsc4py/PETSc/Mat.pyx", line 1525, in petsc4py.PETSc.Mat.ptap petsc4py.PETSc.Error: error code 60 [0] MatPtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9896 [0] MatProductSetFromOptions() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:541 [0] MatProductSetFromOptions_Private() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:435 [0] MatProductSetFromOptions_MPIAIJ() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2372 [0] MatProductSetFromOptions_MPIAIJ_PtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2266 [0] Nonconforming object sizes [0] Matrix local dimensions are incompatible, Acol (0, 100) != Prow (0,34) Abort(1) on node 0 (rank 0 in comm 496): application called MPI_Abort(PYOP2_COMM_WORLD, 1) - process 0 Any thoughts? Thanks, Thanos > On 5 Oct 2023, at 14:23, Thanasis Boutsikakis > <thanasis.boutsika...@corintis.com> wrote: > > This works Pierre. Amazing input, thanks a lot! > >> On 5 Oct 2023, at 14:17, Pierre Jolivet <pie...@joliv.et> wrote: >> >> Not a petsc4py expert here, but you may to try instead: >> A_prime = A.ptap(Phi) >> >> Thanks, >> Pierre >> >>> On 5 Oct 2023, at 2:02 PM, Thanasis Boutsikakis >>> <thanasis.boutsika...@corintis.com> wrote: >>> >>> Thanks Pierre! So I tried this and got a segmentation fault. Is this >>> supposed to work right off the bat or am I missing sth? >>> >>> [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. >>> Abort(59) on node 0 (rank 0 in comm 0): application called >>> MPI_Abort(MPI_COMM_WORLD, 59) - process 0 >>> >>> """Experimenting with PETSc mat-mat multiplication""" >>> >>> import time >>> >>> import numpy as np >>> from colorama import Fore >>> from firedrake import COMM_SELF, COMM_WORLD >>> from firedrake.petsc import PETSc >>> from mpi4py import MPI >>> from numpy.testing import assert_array_almost_equal >>> >>> from utilities import ( >>> Print, >>> create_petsc_matrix, >>> print_matrix_partitioning, >>> ) >>> >>> nproc = COMM_WORLD.size >>> rank = COMM_WORLD.rank >>> >>> # -------------------------------------------- >>> # EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc >>> matrix Phi >>> # A' = Phi.T * A * Phi >>> # [k x k] <- [k x m] x [m x m] x [m x k] >>> # -------------------------------------------- >>> >>> m, k = 11, 7 >>> # 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, m)) >>> Phi_np = np.random.randint(low=0, high=6, size=(m, k)) >>> >>> # -------------------------------------------- >>> # TEST: Galerking projection of numpy matrices A_np and Phi_np >>> # -------------------------------------------- >>> Aprime_np = Phi_np.T @ A_np @ Phi_np >>> Print(f"MATRIX Aprime_np [{Aprime_np.shape[0]}x{Aprime_np.shape[1]}]") >>> Print(f"{Aprime_np}") >>> >>> # Create A as an mpi matrix distributed on each process >>> A = create_petsc_matrix(A_np, sparse=False) >>> >>> # Create Phi as an mpi matrix distributed on each process >>> Phi = create_petsc_matrix(Phi_np, sparse=False) >>> >>> # Create an empty PETSc matrix object to store the result of the PtAP >>> operation. >>> # This will hold the result A' = Phi.T * A * Phi after the computation. >>> A_prime = create_petsc_matrix(np.zeros((k, k)), sparse=False) >>> >>> # Perform the PtAP (Phi Transpose times A times Phi) operation. >>> # In mathematical terms, this operation is A' = Phi.T * A * Phi. >>> # A_prime will store the result of the operation. >>> Phi.PtAP(A, A_prime) >>> >>>> On 5 Oct 2023, at 13:22, Pierre Jolivet <pie...@joliv.et> wrote: >>>> >>>> How about using ptap which will use MatPtAP? >>>> It will be more efficient (and it will help you bypass the issue). >>>> >>>> Thanks, >>>> Pierre >>>> >>>>> On 5 Oct 2023, at 1:18 PM, Thanasis Boutsikakis >>>>> <thanasis.boutsika...@corintis.com> wrote: >>>>> >>>>> Sorry, forgot function create_petsc_matrix() >>>>> >>>>> def create_petsc_matrix(input_array sparse=True): >>>>> """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 >>>>> >>>>> size = ((None, global_rows), (global_cols, global_cols)) >>>>> >>>>> # Create a sparse or dense matrix based on the 'sparse' argument >>>>> if sparse: >>>>> matrix = PETSc.Mat().createAIJ(size=size, comm=COMM_WORLD) >>>>> else: >>>>> matrix = PETSc.Mat().createDense(size=size, comm=COMM_WORLD) >>>>> 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 >>>>> >>>>>> On 5 Oct 2023, at 13:09, Thanasis Boutsikakis >>>>>> <thanasis.boutsika...@corintis.com> wrote: >>>>>> >>>>>> Hi everyone, >>>>>> >>>>>> I am trying a Galerkin projection (see MFE below) and I cannot get the >>>>>> Phi.transposeMatMult(A, A1) work. The error is >>>>>> >>>>>> Phi.transposeMatMult(A, A1) >>>>>> File "petsc4py/PETSc/Mat.pyx", line 1514, in >>>>>> petsc4py.PETSc.Mat.transposeMatMult >>>>>> petsc4py.PETSc.Error: error code 56 >>>>>> [0] MatTransposeMatMult() at >>>>>> /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:10135 >>>>>> [0] MatProduct_Private() at >>>>>> /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9989 >>>>>> [0] No support for this operation for this object type >>>>>> [0] Call MatProductCreate() first >>>>>> >>>>>> Do you know if these exposed to petsc4py or maybe there is another way? >>>>>> I cannot get the MFE to work (neither in sequential nor in parallel) >>>>>> >>>>>> """Experimenting with PETSc mat-mat multiplication""" >>>>>> >>>>>> import time >>>>>> >>>>>> import numpy as np >>>>>> from colorama import Fore >>>>>> from firedrake import COMM_SELF, COMM_WORLD >>>>>> from firedrake.petsc import PETSc >>>>>> from mpi4py import MPI >>>>>> from numpy.testing import assert_array_almost_equal >>>>>> >>>>>> from utilities import ( >>>>>> Print, >>>>>> create_petsc_matrix, >>>>>> ) >>>>>> >>>>>> nproc = COMM_WORLD.size >>>>>> rank = COMM_WORLD.rank >>>>>> >>>>>> # -------------------------------------------- >>>>>> # EXP: Galerkin projection of an mpi PETSc matrix A with an mpi PETSc >>>>>> matrix Phi >>>>>> # A' = Phi.T * A * Phi >>>>>> # [k x k] <- [k x m] x [m x m] x [m x k] >>>>>> # -------------------------------------------- >>>>>> >>>>>> m, k = 11, 7 >>>>>> # 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, m)) >>>>>> Phi_np = np.random.randint(low=0, high=6, size=(m, k)) >>>>>> >>>>>> # Create A as an mpi matrix distributed on each process >>>>>> A = create_petsc_matrix(A_np) >>>>>> >>>>>> # Create Phi as an mpi matrix distributed on each process >>>>>> Phi = create_petsc_matrix(Phi_np) >>>>>> >>>>>> A1 = create_petsc_matrix(np.zeros((k, m))) >>>>>> >>>>>> # Now A1 contains the result of Phi^T * A >>>>>> Phi.transposeMatMult(A, A1) >>>>>> >>>>> >>>> >>> >> >