Pierre, I see your point, but my experiment shows that it does not even run due to size mismatch, so I don’t see how being sparse would change things here. There must be some kind of problem with the parallel ptap(), because it does run sequentially. In order to test that, I changed the flags of the matrix creation to sparse=True and ran it again. Here is the code
"""Experimenting with PETSc mat-mat multiplication""" import numpy as np from firedrake import COMM_WORLD from firedrake.petsc import PETSc 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 # Create A as an mpi matrix distributed on each process A = create_petsc_matrix(A_np, sparse=True) # Create Phi as an mpi matrix distributed on each process Phi = create_petsc_matrix(Phi_np, sparse=True) # 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=True) # 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) I got Traceback (most recent call last): File "/Users/boutsitron/petsc-experiments/mat_vec_multiplication2.py", line 89, in <module> Traceback (most recent call last): File "/Users/boutsitron/petsc-experiments/mat_vec_multiplication2.py", line 89, in <module> Traceback (most recent call last): File "/Users/boutsitron/petsc-experiments/mat_vec_multiplication2.py", line 89, in <module> A_prime = A.ptap(Phi) A_prime = A.ptap(Phi) ^^^^^^^^^^^ File "petsc4py/PETSc/Mat.pyx", line 1525, in petsc4py.PETSc.Mat.ptap A_prime = A.ptap(Phi) ^^^^^^^^^^^ ^^^^^^^^^^^ File "petsc4py/PETSc/Mat.pyx", line 1525, in petsc4py.PETSc.Mat.ptap 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 petsc4py.PETSc.Error: error code 60 [1] MatPtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9896 [1] MatProductSetFromOptions() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:541 [1] MatProductSetFromOptions_Private() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:435 [1] MatProductSetFromOptions_MPIAIJ() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2372 [1] MatProductSetFromOptions_MPIAIJ_PtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2266 [1] Nonconforming object sizes [1] Matrix local dimensions are incompatible, Acol (100, 200) != Prow (34,67) Abort(1) on node 1 (rank 1 in comm 496): application called MPI_Abort(PYOP2_COMM_WORLD, 1) - process 1 petsc4py.PETSc.Error: error code 60 [2] MatPtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matrix.c:9896 [2] MatProductSetFromOptions() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:541 [2] MatProductSetFromOptions_Private() at /Users/boutsitron/firedrake/src/petsc/src/mat/interface/matproduct.c:435 [2] MatProductSetFromOptions_MPIAIJ() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2372 [2] MatProductSetFromOptions_MPIAIJ_PtAP() at /Users/boutsitron/firedrake/src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c:2266 [2] Nonconforming object sizes [2] Matrix local dimensions are incompatible, Acol (200, 300) != Prow (67,100) Abort(1) on node 2 (rank 2 in comm 496): application called MPI_Abort(PYOP2_COMM_WORLD, 1) - process 2 > On 11 Oct 2023, at 07:18, Pierre Jolivet <pie...@joliv.et> wrote: > > I disagree with what Mark and Matt are saying: your code is fine, the error > message is fine, petsc4py is fine (in this instance). > It’s not a typical use case of MatPtAP(), which is mostly designed for > MatAIJ, not MatDense. > On the one hand, in the MatDense case, indeed there will be a mismatch > between the number of columns of A and the number of rows of P, as written in > the error message. > On the other hand, there is not much to optimize when computing C = P’ A P > with everything being dense. > I would just write this as B = A P and then C = P’ B (but then you may face > the same issue as initially reported, please let us know then). > > Thanks, > Pierre > >> On 11 Oct 2023, at 2:42 AM, Mark Adams <mfad...@lbl.gov> wrote: >> >> This looks like a false positive or there is some subtle bug here that we >> are not seeing. >> Could this be the first time parallel PtAP has been used (and reported) in >> petsc4py? >> >> Mark >> >> On Tue, Oct 10, 2023 at 8:27 PM Matthew Knepley <knep...@gmail.com >> <mailto:knep...@gmail.com>> wrote: >>> On Tue, Oct 10, 2023 at 5:34 PM Thanasis Boutsikakis >>> <thanasis.boutsika...@corintis.com >>> <mailto:thanasis.boutsika...@corintis.com>> wrote: >>>> 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). >>> >>> Are you positive about this? P^T A P is designed to run in this scenario, >>> so either we have a bug or the diagnosis is wrong. >>> >>> Thanks, >>> >>> Matt >>> >>>> """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 >>>>> <mailto: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 >>>>>> <mailto: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 >>>>>>> <mailto: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 >>>>>>>> <mailto: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 >>>>>>>>> <mailto: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 >>>>>>>>>> <mailto: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) >>>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>> >>>>>> >>>>> >>>> >>> >>> >>> -- >>> What most experimenters take for granted before they begin their >>> experiments is infinitely more interesting than any results to which their >>> experiments lead. >>> -- Norbert Wiener >>> >>> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/> >