Re: [petsc-users] Orthogonalization of a (sparse) PETSc matrix
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 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
>> 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)
>>
>> #
>> #
Re: [petsc-users] Orthogonalization of a (sparse) PETSc matrix
Ah, there is
https://petsc.org/release/manualpages/Mat/MATSOLVERSPQR/#matsolverspqr See
also https://petsc.org/release/manualpages/Mat/MatGetFactor/#matgetfactor and
https://petsc.org/release/manualpages/Mat/MatQRFactorSymbolic/
> On Aug 29, 2023, at 1:17 PM, Jed Brown wrote:
>
> Suitesparse includes a sparse QR algorithm. The main issue is that (even with
> pivoting) the R factor has the same nonzero structure as a Cholesky factor of
> A^T A, which is generally much denser than a factor of A, and this degraded
> sparsity impacts Q as well.
>
> I wonder if someone would like to contribute a sparse QR to PETSc. It could
> have a default implementation via Cholesky QR and the ability to call SPQR
> from Suitesparse.
>
> Barry Smith writes:
>
>> Are the nonzero structures of all the rows related? If they are, one could
>> devise a routine to take advantage of this relationship, but if the nonzero
>> structures of each row are "randomly" different from all the other rows,
>> then it is difficult to see how one can take advantage of the sparsity.
>>
>>
>>
>>> On Aug 29, 2023, at 12:50 PM, Thanasis Boutsikakis
>>> wrote:
>>>
>>> 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
>>>
Re: [petsc-users] Error while building PETSc with MATLAB
Well - you sent in libmesh log not petsc's configure.log/make.log for petsc-3.17 Anyway - with petsc-3.13 - you have: Matlab: Includes: -I/usr/local/MATLAB/R2020b/extern/include /usr/local/MATLAB/R2020b MatlabEngine: Library: -Wl,-rpath,/usr/local/MATLAB/R2020b/sys/os/glnxa64:/usr/local/MATLAB/R2020b/bin/glnxa64:/usr/local/MATLAB/R2020b/extern/lib/glnxa64 -L/usr/local/MATLAB/R2020b/bin/glnxa64 -L/usr/local/MATLAB/R2020b/extern/lib/glnxa64 -leng -lmex -lmx -lmat -lut -licudata -licui18n -licuuc Language used to compile PETSc: C < With petsc-3.19 (and matlab-R2022a) - we are seeing: https://gitlab.com/petsc/petsc/-/jobs/4904566768 >>> Matlab: Includes: -I/nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a/extern/include Libraries: -Wl,-rpath,/nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a/bin/glnxa64 -L/nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a/bin/glnxa64 -leng -lmex -lmx -lmat Executable: /nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a mex: /nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a/bin/mex matlab: /nfs/gce/software/custom/linux-ubuntu22.04-x86_64/matlab/R2022a/bin/matlab -glnxa64 <<< I.e "-lut -licudata -licui18n -licuuc" are not preset here. This might be a change wrt newer matlab versions. You can: - edit /home/vit/sfw/petsc/3.13.4/linux-opt/lib/petsc/conf/petscvariables and remove all occurrences of "-lut -licudata -licui18n -licuuc" - now run 'make all' in '/home/vit/sfw/petsc/3.13.4' And see if the build works now. Satish On Tue, 29 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > I am sharing the log files while building petsc3.13.4 with matlab and also > the log file while building libmesh with petsc3.17.5 and matlab. Building > petsc 3.17.5 with matlab was done successfully. But libmesh is not able to > find the petsc > Please find the attachments. > > On Tue, Aug 29, 2023 at 7:31 PM Satish Balay wrote: > > > Send configure.log, make.log from both petsc-3.13 and 3.17 [or 3.19]. > > > > [you can gzip them to make the logs friendly to mailing list - or send > > them to petsc-maint] > > > > And does test suite work with 3.17? [or 3.19?] > > > > Satish > > > > On Tue, 29 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > > > > > I am sharing the make.log file while building petsc-3.13.4 with Matlab. > > > Please find the attachment and do the needful. > > > > > > On Tue, Aug 29, 2023 at 10:19 AM INTURU SRINIVAS 20PHD0548 < > > > inturu.srinivas2...@vitstudent.ac.in> wrote: > > > > > > > I tried with petsc-3.17.5. During building of libmesh, the error shows > > > > petsc was not found > > > > > > > > On Mon, Aug 28, 2023 at 9:43 PM Satish Balay > > wrote: > > > > > > > >> https://ibamr.github.io/linux says petsc-3.17 > > > >> > > > >> Here you are using 3.13 > > > >> > > > >> Can you retry with petsc-3.17.5? > > > >> > > > >> Satish > > > >> > > > >> On Mon, 28 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > > > >> > > > >> > Hello, > > > >> > > > > >> > I want to build PETSc with MATLAB for working on the simulation > > using > > > >> IBAMR > > > >> > open software. While building the PETSc, using the following > > > >> > > > > >> > export PETSC_DIR=$PWD > > > >> > export PETSC_ARCH=linux-debug > > > >> > ./configure \ > > > >> > --CC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicc \ > > > >> > --CXX=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicxx \ > > > >> > --FC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpif90 \ > > > >> > --with-debugging=1 \ > > > >> > --download-hypre=1 \ > > > >> > --download-fblaslapack=1 \ > > > >> > --with-x=0 \ > > > >> > --with-matlab-dir=/usr/local/MATLAB/R2020b/ > > > >> > --with-matlab-engine=1 > > > >> > --with-matlab-engine-dir=/usr/local/MATLAB/R2020b/extern/engines/ > > > >> > > > > >> > make -j4 > > > >> > make -j4 test > > > >> > > > > >> > I got the following error > > > >> > CLINKER > > > >> linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test > > > >> > /usr/bin/ld: > > > >> > > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > > >> > function `EvaluateResidual': > > > >> > > > > >> > > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:32: > > > >> > undefined reference to `PetscMatlabEnginePut' > > > >> > /usr/bin/ld: > > > >> > > > > >> > > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:33: > > > >> > undefined reference to `PetscMatlabEngineEvaluate' > > > >> > /usr/bin/ld: > > > >> > > > > >> > > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:35: > > > >> > undefined reference to `PetscMatlabEngineGet' > > > >> > /usr/bin/ld: > > > >> > > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > > >> > function `EvaluateJacobian': > > > >> > > > > >> > >
Re: [petsc-users] Orthogonalization of a (sparse) PETSc matrix
Suitesparse includes a sparse QR algorithm. The main issue is that (even with
pivoting) the R factor has the same nonzero structure as a Cholesky factor of
A^T A, which is generally much denser than a factor of A, and this degraded
sparsity impacts Q as well.
I wonder if someone would like to contribute a sparse QR to PETSc. It could
have a default implementation via Cholesky QR and the ability to call SPQR from
Suitesparse.
Barry Smith writes:
> Are the nonzero structures of all the rows related? If they are, one could
> devise a routine to take advantage of this relationship, but if the nonzero
> structures of each row are "randomly" different from all the other rows, then
> it is difficult to see how one can take advantage of the sparsity.
>
>
>
>> On Aug 29, 2023, at 12:50 PM, Thanasis Boutsikakis
>> wrote:
>>
>> 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
>>
Re: [petsc-users] Orthogonalization of a (sparse) PETSc matrix
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
> 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
Re: [petsc-users] Orthogonalization of a (sparse) PETSc matrix
Are the nonzero structures of all the rows related? If they are, one could
devise a routine to take advantage of this relationship, but if the nonzero
structures of each row are "randomly" different from all the other rows, then
it is difficult to see how one can take advantage of the sparsity.
> On Aug 29, 2023, at 12:50 PM, Thanasis Boutsikakis
> wrote:
>
> 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),
Re: [petsc-users] Error while building PETSc with MATLAB
On Tue, Aug 29, 2023 at 9:08 AM Satish Balay via petsc-users < petsc-users@mcs.anl.gov> wrote: > Send configure.log, make.log from both petsc-3.13 and 3.17 [or 3.19]. > > [you can gzip them to make the logs friendly to mailing list - or send > them to petsc-maint] > > And does test suite work with 3.17? [or 3.19?] > David Wells is working on this. The change is that petscversion.h now includes petscconf.h which means you need all the include flags, but Libmesh does not get the flags right. Thanks, Matt > Satish > > On Tue, 29 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > > > I am sharing the make.log file while building petsc-3.13.4 with Matlab. > > Please find the attachment and do the needful. > > > > On Tue, Aug 29, 2023 at 10:19 AM INTURU SRINIVAS 20PHD0548 < > > inturu.srinivas2...@vitstudent.ac.in> wrote: > > > > > I tried with petsc-3.17.5. During building of libmesh, the error shows > > > petsc was not found > > > > > > On Mon, Aug 28, 2023 at 9:43 PM Satish Balay > wrote: > > > > > >> https://ibamr.github.io/linux says petsc-3.17 > > >> > > >> Here you are using 3.13 > > >> > > >> Can you retry with petsc-3.17.5? > > >> > > >> Satish > > >> > > >> On Mon, 28 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > > >> > > >> > Hello, > > >> > > > >> > I want to build PETSc with MATLAB for working on the simulation > using > > >> IBAMR > > >> > open software. While building the PETSc, using the following > > >> > > > >> > export PETSC_DIR=$PWD > > >> > export PETSC_ARCH=linux-debug > > >> > ./configure \ > > >> > --CC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicc \ > > >> > --CXX=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicxx \ > > >> > --FC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpif90 \ > > >> > --with-debugging=1 \ > > >> > --download-hypre=1 \ > > >> > --download-fblaslapack=1 \ > > >> > --with-x=0 \ > > >> > --with-matlab-dir=/usr/local/MATLAB/R2020b/ > > >> > --with-matlab-engine=1 > > >> > --with-matlab-engine-dir=/usr/local/MATLAB/R2020b/extern/engines/ > > >> > > > >> > make -j4 > > >> > make -j4 test > > >> > > > >> > I got the following error > > >> > CLINKER > > >> linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test > > >> > /usr/bin/ld: > > >> > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > >> > function `EvaluateResidual': > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:32: > > >> > undefined reference to `PetscMatlabEnginePut' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:33: > > >> > undefined reference to `PetscMatlabEngineEvaluate' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:35: > > >> > undefined reference to `PetscMatlabEngineGet' > > >> > /usr/bin/ld: > > >> > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > >> > function `EvaluateJacobian': > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:46: > > >> > undefined reference to `PetscMatlabEnginePut' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:47: > > >> > undefined reference to `PetscMatlabEngineEvaluate' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:49: > > >> > undefined reference to `PetscMatlabEngineGet' > > >> > /usr/bin/ld: > > >> > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > >> > function `TaoPounders': > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:75: > > >> > undefined reference to `PetscMatlabEngineGet' > > >> > /usr/bin/ld: > > >> > > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > > >> > function `main': > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:126: > > >> > undefined reference to `PetscMatlabEngineCreate' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:127: > > >> > undefined reference to `PetscMatlabEngineEvaluate' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:139: > > >> > undefined reference to `PetscMatlabEngineEvaluate' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:140: > > >> > undefined reference to `PetscMatlabEngineGetArray' > > >> > /usr/bin/ld: > > >> > > > >> > /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:142: > > >> > undefined reference to `PetscMatlabEngineGetArray' > > >> > /usr/bin/ld: > > >> > > > >> >
[petsc-users] Orthogonalization of a (sparse) PETSc matrix
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,
)
Re: [petsc-users] Error while building PETSc with MATLAB
Send configure.log, make.log from both petsc-3.13 and 3.17 [or 3.19]. [you can gzip them to make the logs friendly to mailing list - or send them to petsc-maint] And does test suite work with 3.17? [or 3.19?] Satish On Tue, 29 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > I am sharing the make.log file while building petsc-3.13.4 with Matlab. > Please find the attachment and do the needful. > > On Tue, Aug 29, 2023 at 10:19 AM INTURU SRINIVAS 20PHD0548 < > inturu.srinivas2...@vitstudent.ac.in> wrote: > > > I tried with petsc-3.17.5. During building of libmesh, the error shows > > petsc was not found > > > > On Mon, Aug 28, 2023 at 9:43 PM Satish Balay wrote: > > > >> https://ibamr.github.io/linux says petsc-3.17 > >> > >> Here you are using 3.13 > >> > >> Can you retry with petsc-3.17.5? > >> > >> Satish > >> > >> On Mon, 28 Aug 2023, INTURU SRINIVAS 20PHD0548 via petsc-users wrote: > >> > >> > Hello, > >> > > >> > I want to build PETSc with MATLAB for working on the simulation using > >> IBAMR > >> > open software. While building the PETSc, using the following > >> > > >> > export PETSC_DIR=$PWD > >> > export PETSC_ARCH=linux-debug > >> > ./configure \ > >> > --CC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicc \ > >> > --CXX=$HOME/sfw/linux/openmpi/4.1.4/bin/mpicxx \ > >> > --FC=$HOME/sfw/linux/openmpi/4.1.4/bin/mpif90 \ > >> > --with-debugging=1 \ > >> > --download-hypre=1 \ > >> > --download-fblaslapack=1 \ > >> > --with-x=0 \ > >> > --with-matlab-dir=/usr/local/MATLAB/R2020b/ > >> > --with-matlab-engine=1 > >> > --with-matlab-engine-dir=/usr/local/MATLAB/R2020b/extern/engines/ > >> > > >> > make -j4 > >> > make -j4 test > >> > > >> > I got the following error > >> > CLINKER > >> linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test > >> > /usr/bin/ld: > >> > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > >> > function `EvaluateResidual': > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:32: > >> > undefined reference to `PetscMatlabEnginePut' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:33: > >> > undefined reference to `PetscMatlabEngineEvaluate' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:35: > >> > undefined reference to `PetscMatlabEngineGet' > >> > /usr/bin/ld: > >> > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > >> > function `EvaluateJacobian': > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:46: > >> > undefined reference to `PetscMatlabEnginePut' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:47: > >> > undefined reference to `PetscMatlabEngineEvaluate' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:49: > >> > undefined reference to `PetscMatlabEngineGet' > >> > /usr/bin/ld: > >> > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > >> > function `TaoPounders': > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:75: > >> > undefined reference to `PetscMatlabEngineGet' > >> > /usr/bin/ld: > >> > linux-debug/tests/tao/leastsquares/tutorials/matlab/matlab_ls_test.o: in > >> > function `main': > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:126: > >> > undefined reference to `PetscMatlabEngineCreate' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:127: > >> > undefined reference to `PetscMatlabEngineEvaluate' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:139: > >> > undefined reference to `PetscMatlabEngineEvaluate' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:140: > >> > undefined reference to `PetscMatlabEngineGetArray' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:142: > >> > undefined reference to `PetscMatlabEngineGetArray' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:144: > >> > undefined reference to `PetscMatlabEngineGetArray' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:146: > >> > undefined reference to `PetscMatlabEngineGetArray' > >> > /usr/bin/ld: > >> > > >> /home/vit/sfw/petsc/3.13.4/src/tao/leastsquares/tutorials/matlab/matlab_ls_test.c:148: > >> > undefined reference to `PetscMatlabEngineGetArray' > >> > /usr/bin/ld: > >> > > >>
