> On 11 Oct 2023, at 9:13 AM, Thanasis Boutsikakis > <[email protected]> wrote: > > 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?
Hopefully this should make things clearer to you: https://petsc.org/release/manual/mat/#sec-matlayout Thanks, Pierre > 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/> >>>> >>> >> >
