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) >>>> >>> >> >