Very good catch Pierre, thanks a lot! This made everything work: the two-step process and the ptap(). I mistakenly thought that I should not let the local number of columns to be None, since the matrix is only partitioned row-wise. Could you please explain what happened because of my setting the local column number so that I get the philosophy behind this partitioning?
Thanks again, Thanos > On 11 Oct 2023, at 09:04, Pierre Jolivet <[email protected]> wrote: > > That’s because: > size = ((None, global_rows), (global_cols, global_cols)) > should be: > size = ((None, global_rows), (None, global_cols)) > Then, it will work. > $ ~/repo/petsc/arch-darwin-c-debug-real/bin/mpirun -n 4 python3.12 test.py && > echo $? > 0 > > Thanks, > Pierre > >> On 11 Oct 2023, at 8:58 AM, Thanasis Boutsikakis >> <[email protected]> wrote: >> >> 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 <[email protected]> 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 <[email protected]> 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 <[email protected] >>>> <mailto:[email protected]>> wrote: >>>>> On Tue, Oct 10, 2023 at 5:34 PM Thanasis Boutsikakis >>>>> <[email protected] >>>>> <mailto:[email protected]>> 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 >>>>>>> <[email protected] >>>>>>> <mailto:[email protected]>> wrote: >>>>>>> >>>>>>> This works Pierre. Amazing input, thanks a lot! >>>>>>> >>>>>>>> On 5 Oct 2023, at 14:17, Pierre Jolivet <[email protected] >>>>>>>> <mailto:[email protected]>> 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 >>>>>>>>> <[email protected] >>>>>>>>> <mailto:[email protected]>> 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 <[email protected] >>>>>>>>>> <mailto:[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] >>>>>>>>>>> <mailto:[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] >>>>>>>>>>>> <mailto:[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) >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>> >>>>>> >>>>> >>>>> >>>>> -- >>>>> 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/> >>> >> >
