[petsc-users] Galerkin projection using petsc4py
Pierre Jolivet
pierre at joliv.et
Wed Oct 11 02:04:51 CDT 2023
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/17150218/attachment-0001.html>
More information about the petsc-users
mailing list