[petsc-users] Galerkin projection using petsc4py

Thanasis Boutsikakis thanasis.boutsikakis at corintis.com
Wed Oct 11 01:58:18 CDT 2023 

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/c80d6864/attachment-0001.html>

More information about the petsc-users mailing list