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