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