[petsc-users] Galerkin projection using petsc4py

Thanasis Boutsikakis thanasis.boutsikakis at corintis.com
Wed Oct 11 02:13:28 CDT 2023 

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 <pierre at joliv.et> 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 <thanasis.boutsikakis at corintis.com> 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 <pierre at 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 <mfadams at 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 <knepley at gmail.com <mailto:knepley at gmail.com>> wrote:
>>>>> On Tue, Oct 10, 2023 at 5:34 PM Thanasis Boutsikakis <thanasis.boutsikakis at corintis.com <mailto:thanasis.boutsikakis at 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.boutsikakis at corintis.com <mailto:thanasis.boutsikakis at corintis.com>> wrote:
>>>>>>> 
>>>>>>> This works Pierre. Amazing input, thanks a lot!
>>>>>>> 
>>>>>>>> On 5 Oct 2023, at 14:17, Pierre Jolivet <pierre at joliv.et <mailto:pierre at 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.boutsikakis at corintis.com <mailto:thanasis.boutsikakis at 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 <pierre at joliv.et <mailto:pierre at 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.boutsikakis at corintis.com <mailto:thanasis.boutsikakis at 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.boutsikakis at corintis.com <mailto:thanasis.boutsikakis at 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/>
>>> 
>> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231011/b85869ef/attachment-0001.html>

More information about the petsc-users mailing list