[petsc-users] Galerkin projection using petsc4py

Thu Oct 5 07:17:52 CDT 2023

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> 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> 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> 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> 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)
>>>> 
>>> 
>> 
> 

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