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

[Jose E. Roman]

The conversion from MATAIJ to MATDENSE should be very cheap, see https://gitlab.com/petsc/petsc/-/blob/main/src/mat/impls/dense/seq/dense.c?ref_type=heads#L172

The matrix copy hidden inside createFromMat() is likely more expensive. I am currently working on a modification of BV that will be included in version 3.20 if everything goes well - then I think I can allow passing a sparse matrix to createFromMat() and do the conversion internally, avoiding the matrix copy.

Jose


> El 29 ago 2023, a las 22:46, Thanasis Boutsikakis <thanasis.boutsikakis at corintis.com> escribió:
> 
> 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,
>>> )
>> 
>