[petsc-users] Orthogonalization of a (sparse) PETSc matrix

[Thanasis Boutsikakis]

Thanks Jose, 

This works indeed. However, I was under the impression that this conversion might be very costly for big matrices with low sparsity and it would scale with the number of non-zero values.

Do you have any idea of the efficiency of this operation?

Thanks

> On 29 Aug 2023, at 19:13, Jose E. Roman <jroman at dsic.upv.es> wrote:
> 
> The result of bv.orthogonalize() is most probably a dense matrix, and the result replaces the input matrix, that's why the input matrix is required to be dense.
> 
> You can simply do this:
> 
>  bv = SLEPc.BV().createFromMat(A.convert('dense'))
> 
> Jose
> 
>> El 29 ago 2023, a las 18:50, Thanasis Boutsikakis <thanasis.boutsikakis at corintis.com> escribió:
>> 
>> Hi all, I have the following code that orthogonalizes a PETSc matrix. The problem is that this implementation requires that the PETSc matrix is dense, otherwise, it fails at bv.SetFromOptions(). Hence the assert in orthogonality().
>> 
>> What could I do in order to be able to orthogonalize sparse matrices as well? Could I convert it efficiently? (I tried to no avail)
>> 
>> Thanks!
>> 
>> """Experimenting with matrix orthogonalization"""
>> 
>> import contextlib
>> import sys
>> import time
>> import numpy as np
>> from firedrake import COMM_WORLD
>> from firedrake.petsc import PETSc
>> 
>> import slepc4py
>> 
>> slepc4py.init(sys.argv)
>> from slepc4py import SLEPc
>> 
>> from numpy.testing import assert_array_almost_equal
>> 
>> EPSILON_USER = 1e-4
>> EPS = sys.float_info.epsilon
>> 
>> 
>> def Print(message: str):
>>    """Print function that prints only on rank 0 with color
>> 
>>    Args:
>>        message (str): message to be printed
>>    """
>>    PETSc.Sys.Print(message)
>> 
>> 
>> 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
>> 
>> 
>> def orthogonality(A):  # sourcery skip: avoid-builtin-shadow
>>    """Checking and correcting orthogonality
>> 
>>    Args:
>>        A (PETSc.Mat): Matrix of size [m x k].
>> 
>>    Returns:
>>        PETSc.Mat: Matrix of size [m x k].
>>    """
>>    # Check if the matrix is dense
>>    mat_type = A.getType()
>>    assert mat_type in (
>>        "seqdense",
>>        "mpidense",
>>    ), "A must be a dense matrix. SLEPc.BV().createFromMat() requires a dense matrix."
>> 
>>    m, k = A.getSize()
>> 
>>    Phi1 = A.getColumnVector(0)
>>    Phi2 = A.getColumnVector(k - 1)
>> 
>>    # Compute dot product using PETSc function
>>    dot_product = Phi1.dot(Phi2)
>> 
>>    if abs(dot_product) > min(EPSILON_USER, EPS * m):
>>        Print("    Matrix is not orthogonal")
>> 
>>        # Type can be CHOL, GS, mro(), SVQB, TSQR, TSQRCHOL
>>        _type = SLEPc.BV().OrthogBlockType.GS
>> 
>>        bv = SLEPc.BV().createFromMat(A)
>>        bv.setFromOptions()
>>        bv.setOrthogonalization(_type)
>>        bv.orthogonalize()
>> 
>>        A = bv.createMat()
>> 
>>        Print("    Matrix successfully orthogonalized")
>> 
>>        # # Assembly the matrix to compute the final structure
>>        if not A.assembled:
>>            A.assemblyBegin()
>>            A.assemblyEnd()
>>    else:
>>        Print("    Matrix is orthogonal")
>> 
>>    return A
>> 
>> 
>> # --------------------------------------------
>> # EXP: Orthogonalization of an mpi PETSc matrix
>> # --------------------------------------------
>> 
>> 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, k))
>> 
>> A = create_petsc_matrix(A_np, sparse=False)
>> 
>> A_orthogonal = orthogonality(A)
>> 
>> # --------------------------------------------
>> # TEST: Orthogonalization of a numpy matrix
>> # --------------------------------------------
>> # Generate A_np_orthogonal
>> A_np_orthogonal, _ = np.linalg.qr(A_np)
>> 
>> # Get the local values from A_orthogonal
>> local_rows_start, local_rows_end = A_orthogonal.getOwnershipRange()
>> A_orthogonal_local = A_orthogonal.getValues(
>>    range(local_rows_start, local_rows_end), range(k)
>> )
>> 
>> # Assert the correctness of the multiplication for the local subset
>> assert_array_almost_equal(
>>    np.abs(A_orthogonal_local),
>>    np.abs(A_np_orthogonal[local_rows_start:local_rows_end, :]),
>>    decimal=5,
>> )
>