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 <[email protected]> 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 >> <[email protected]> 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 >>> <[email protected]> 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) >>> >> >
