January 2024 - petsc-users - mailman3.cels.anl.gov

ScaLAPACK EPS error
by Peder Jørgensgaard Olesen 19 Jan '24

19 Jan '24

Hello, I need to determine the full set of eigenpairs to a rather large (N=16,000) dense Hermitian matrix. I've managed to do this using SLEPc's standard Krylov-Schur EPS, but I think it could be done more efficiently using ScaLAPACK. I receive the following error when attempting this. As I understand it, descinit is used to initialize an array, and the variable in question designates the leading dimension of the array, for which it seems an illegal value is somehow passed. I know ScaLAPACK is an external package, but it seems as if the error would be in the call from SLEPc. Any ideas as to what could cause this? Thanks, Peder Error message (excerpt): PETSC ERROR: #1 MatConvert_Dense_ScaLAPACK() at [...]/matscalapack.c:1032 PETSC ERROR: #2 MatConvert at [...]/matrix.c:4250 PETSC ERROR: #3 EPSSetUp_ScaLAPACK() at [...]/scalapack.c:47 PETSC ERROR: #4 EPSSetUp() at [...]/epssetup.c:323 PETSC ERROR: #5 EPSSolve at [...]/epssolve.c:134 PETSC ERROR: ------ Error message ------ PETSC ERROR: Error in external library PETSC ERROR: Error in ScaLAPACK subroutine descinit: info=-9 (...) Log file (excerpt): { 357, 0}: On entry to DESCINIT parameter number 9 had an illegal value [and a few hundred lines similar to this]

3 6

MatAssemblyBegin freezes during MPI communication
by 袁煕 18 Jan '24

18 Jan '24

Dear PETSc Experts, My FEM program works well generally, but in some specific cases with multiple CPUs are used, it freezes when calling MatAssemblyBegin where PMPI_Allreduce is called (see attached file). After some investigation, I found that it is most probably due to ・ MatSetValue is not called from all CPUs before MatAssemblyBegin For example, when 4 CPUs are used, if there are elements in CPU 0,1,2 but no elements in CPU 3, then all CPUs other than CPU 3 would call MatSetValue function. I want to know 1. If my conjecture could be right? And If so 2. Are there any convenient means to avoid this problem? Thanks, Xi YUAN, PhD Solid Mechanics

3 3

HashMap Error when populating AIJCUSPARSE matrix
by Yesypenko, Anna 18 Jan '24

18 Jan '24

Dear Petsc users/developers, I'm experiencing a bug when using petsc4py with GPU support. It may be my mistake in how I set up a AIJCUSPARSE matrix. For larger matrices, I sometimes encounter a error in assigning matrix values; the error is thrown in PetscHMapIJVQuerySet(). Here is a minimum snippet that populates a sparse tridiagonal matrix. ``` from petsc4py import PETSc from scipy.sparse import diags import numpy as np n = int(5e5); nnz = 3 * np.ones(n, dtype=np.int32); nnz[0] = nnz[-1] = 2 A = PETSc.Mat(comm=PETSc.COMM_WORLD) A.createAIJ(size=[n,n],comm=PETSc.COMM_WORLD,nnz=nnz) A.setType('aijcusparse') tmp = diags([-1,2,-1],[-1,0,+1],shape=(n,n)).tocsr() A.setValuesCSR(tmp.indptr,tmp.indices,tmp.data) ####### this is the line where the error is thrown. A.assemble() ``` The error trace is below: ``` File "petsc4py/PETSc/Mat.pyx", line 2603, in petsc4py.PETSc.Mat.setValuesCSR File "petsc4py/PETSc/petscmat.pxi", line 1039, in petsc4py.PETSc.matsetvalues_csr File "petsc4py/PETSc/petscmat.pxi", line 1032, in petsc4py.PETSc.matsetvalues_ijv petsc4py.PETSc.Error: error code 76 [0] MatSetValues() at /work/06368/annayesy/ls6/petsc/src/mat/interface/matrix.c:1497 [0] MatSetValues_Seq_Hash() at /work/06368/annayesy/ls6/petsc/include/../src/mat/impls/aij/seq/seqhashmatsetvalues.h:52 [0] PetscHMapIJVQuerySet() at /work/06368/annayesy/ls6/petsc/include/petsc/private/hashmapijv.h:10 [0] Error in external library [0] [khash] Assertion: `ret >= 0' failed. ``` If I run the same script a handful of times, it will run without errors eventually. Does anyone have insight on why it is behaving this way? I'm running on a node with 3x NVIDIA A100 PCIE 40GB. Thank you! Anna

3 10

undefined reference to `petsc_allreduce_ct_th'
by Aaron Scheinberg 18 Jan '24

18 Jan '24

Hello, I'm getting this error when linking: undefined reference to `petsc_allreduce_ct_th' The instances are regular MPI_Allreduces in my code that are not located in parts of the code related to PETSc, so I'm wondering what is happening to involve PETSc here? Can I configure it to avoid that? I consulted google, the FAQ and skimmed other documentation but didn't see anything. Thanks! Aaron

3 3

Undestanding how to increase the overlap
by TARDIEU Nicolas 18 Jan '24

18 Jan '24

Dear PETSc Team, I am trying to understand how to increase the overlap of a matrix. I wrote the attached petsc4py script where I build a simple matrix and play with the increaseOverlap method. Unfortunately, before and after the call, nothing changes in the index set. FYI, I have tried to mimic src/ksp/ksp/tutorials/ex82.c<https://petsc.org/main/src/ksp/ksp/tutorials/ex82.c.html> line 72:76. Here is how I run the script : "mpiexec -n 2 python test_overlap.py" Could you please indicate what I am missing ? Regards, Nicolas -- Nicolas Tardieu Ing PhD Computational Mechanics EDF - R&D Dpt ERMES PARIS-SACLAY, FRANCE Ce message et toutes les pièces jointes (ci-après le 'Message') sont établis à l'intention exclusive des destinataires et les informations qui y figurent sont strictement confidentielles. Toute utilisation de ce Message non conforme à sa destination, toute diffusion ou toute publication totale ou partielle, est interdite sauf autorisation expresse. Si vous n'êtes pas le destinataire de ce Message, il vous est interdit de le copier, de le faire suivre, de le divulguer ou d'en utiliser tout ou partie. Si vous avez reçu ce Message par erreur, merci de le supprimer de votre système, ainsi que toutes ses copies, et de n'en garder aucune trace sur quelque support que ce soit. Nous vous remercions également d'en avertir immédiatement l'expéditeur par retour du message. Il est impossible de garantir que les communications par messagerie électronique arrivent en temps utile, sont sécurisées ou dénuées de toute erreur ou virus. ____________________________________________________ This message and any attachments (the 'Message') are intended solely for the addressees. The information contained in this Message is confidential. Any use of information contained in this Message not in accord with its purpose, any dissemination or disclosure, either whole or partial, is prohibited except formal approval. If you are not the addressee, you may not copy, forward, disclose or use any part of it. If you have received this message in error, please delete it and all copies from your system and notify the sender immediately by return message. E-mail communication cannot be guaranteed to be timely secure, error or virus-free.

2 5

Re: [petsc-users] SNES seems not use my matrix-free operation
by Barry Smith 16 Jan '24

16 Jan '24

> However, my GK_op (which is the reloaded MATOP_MULT) gives me some problem. It is entered but then crashed with Segmentation Violation error. So I guess my indices may be wrong. I wonder do I need to use the local Vec (of dF), and should my output Vec also be in the correct shape (i.e. after calculation I need to transform back into a Vec)? As you can see here, my dF is a tensor defined on every grid point. > The input for the matrix-vector product is a global vector, as is the result. (Not like the arguments to DMSNESSetJacobianLocal). This means that your MATOP_MULT function needs to do the DMGlobalToLocal() vector operation first then the "unwrapping" from the vector to the array format at the beginning of the routine. Similarly it needs to "unwrap" the result vector as an array. See src/snes/tutorials/ex14f.F90 and in particular the code block PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)) PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)) ! Get pointers to vector data PetscCall(VecGetArrayReadF90(localX,xx,ierr)) PetscCall(VecGetArrayF90(F,ff,ierr)) Barry You really shouldn't be using DMSNESSetJacobianLocal() for your code. Basically all the DMSNESSetJacobianLocal() gives you is that it automatically handles the global to local mapping and unwrapping of the vector to an array, but it doesn't work for shell matrices. > On Jan 9, 2024, at 6:30 AM, Yi Hu <y.hu(a)mpie.de> wrote: > > Dear Barry, > > Thanks for your help. > > It works when doing first SNESSetJacobian() with my created shell matrix Jac in the main (or module) and then DMSNESSetJacobianLocal() to associate with my DM and an dummy formJacobian callback (which is doing nothing). My SNES can now recognize my shell matrix and do my customized operation. > > However, my GK_op (which is the reloaded MATOP_MULT) gives me some problem. It is entered but then crashed with Segmentation Violation error. So I guess my indices may be wrong. I wonder do I need to use the local Vec (of dF), and should my output Vec also be in the correct shape (i.e. after calculation I need to transform back into a Vec)? As you can see here, my dF is a tensor defined on every grid point. > > Best wishes, > Yi > > From: Barry Smith <bsmith(a)petsc.dev <mailto:[email protected]>> > Sent: Monday, January 8, 2024 6:41 PM > To: Yi Hu <y.hu(a)mpie.de <mailto:[email protected]>> > Cc: petsc-users(a)mcs.anl.gov <mailto:[email protected]> > Subject: Re: [petsc-users] SNES seems not use my matrix-free operation > > > "formJacobian" should not be __creating__ the matrices. Here "form" means computing the numerical values in the matrix (or when using a shell matrix it means keeping a copy of X so that your custom matrix-free multiply knows the base location where the matrix free Jacobian-vector products are computed.) > > You create the shell matrices up in your main program and pass them in with SNESSetJacobian(). > > Try first calling SNESSetJacobian() to provide the matrices (provide a dummy function argument) and then call DMSNESSetJacobianLocal() to provide your "formjacobian" function (that does not create the matrices). > > Barry > > > Yes, "form" is a bad word that should not have been used in our code. > > > > > On Jan 8, 2024, at 12:24 PM, Yi Hu <y.hu(a)mpie.de <mailto:[email protected]>> wrote: > > Dear PETSc Experts, > > I am implementing a matrix-free jacobian for my SNES solver in Fortran. (command line option -snes_type newtonls -ksp_type gmres) > > In the main program, I define my residual and jacobian and matrix-free jacobian like the following, > > … > call DMDASNESSetFunctionLocal(DM_mech, INSERT_VALUES, formResidual, PETSC_NULL_SNES, err_PETSc) > call DMSNESSetJacobianLocal(DM_mech, formJacobian, PETSC_NULL_SNES, err_PETSc) > … > > subroutine formJacobian(residual_subdomain,F,Jac_pre,Jac,dummy,err_PETSc) > > #include <petsc/finclude/petscmat.h> > use petscmat > implicit None > DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & > residual_subdomain !< DMDA info (needs to be named "in" for macros like XRANGE to work) > real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & > F !< deformation gradient field > Mat :: Jac, Jac_pre > PetscObject :: dummy > PetscErrorCode :: err_PETSc > PetscInt :: N_dof ! global number of DoF, maybe only a placeholder > > N_dof = 9*product(cells(1:2))*cells3 > > print*, 'in my jac' > > call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac,err_PETSc) > CHKERRQ(err_PETSc) > call MatShellSetOperation(Jac,MATOP_MULT,GK_op,err_PETSc) > CHKERRQ(err_PETSc) > > print*, 'in my jac' > > ! for jac preconditioner > call MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N_dof,N_dof,0,Jac_pre,err_PETSc) > CHKERRQ(err_PETSc) > call MatShellSetOperation(Jac_pre,MATOP_MULT,GK_op,err_PETSc) > CHKERRQ(err_PETSc) > > print*, 'in my jac' > > end subroutine formJacobian > > subroutine GK_op(Jac,dF,output,err_PETSc) > real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: & > dF > real(pREAL), dimension(3,3,cells(1),cells(2),cells3), intent(out) :: & > output > real(pREAL), dimension(3,3) :: & > deltaF_aim = 0.0_pREAL > > Mat :: Jac > PetscErrorCode :: err_PETSc > > integer :: i, j, k, e > > … a lot of calculations … > > print*, 'in GK op' > > end subroutine GK_op > > The first question is that: it seems I still need to explicitly define the interface of MatCreateShell() and MatShellSetOperation() to properly use them, even though I include them via “use petscmat”. It is a little bit strange to me, since some examples do not perform this step. > > Then the main issue is that I can build my own Jacobian from my call back function formJacobian, and confirm my Jacobian is a shell matrix (by MatView). However, my customized operator GK_op is not called when solving the nonlinear system (not print my “in GK op”). When I try to monitor my SNES, it gives me some conventional output not mentioning my matrix-free operations. So I guess my customized MATOP_MULT may be not associated with Jacobian. Or my configuration is somehow wrong. Could you help me solve this issue? > > Thanks, > Yi > > > > ------------------------------------------------- > Stay up to date and follow us on LinkedIn, Twitter and YouTube. > > Max-Planck-Institut für Eisenforschung GmbH > Max-Planck-Straße 1 > D-40237 Düsseldorf > > Handelsregister B 2533 > Amtsgericht Düsseldorf > > Geschäftsführung > Prof. Dr. Gerhard Dehm > Prof. Dr. Jörg Neugebauer > Prof. Dr. Dierk Raabe > Dr. Kai de Weldige > > Ust.-Id.-Nr.: DE 11 93 58 514 > Steuernummer: 105 5891 1000 > > > Please consider that invitations and e-mails of our institute are > only valid if they end with …(a)mpie.de. > If you are not sure of the validity please contact rco(a)mpie.de <mailto:[email protected]> > > Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails > aus unserem Haus nur mit der Endung …(a)mpie.de gültig sind. > In Zweifelsfällen wenden Sie sich bitte an rco(a)mpie.de <mailto:[email protected]> > ------------------------------------------------- > > > > ------------------------------------------------- > Stay up to date and follow us on LinkedIn, Twitter and YouTube. > > Max-Planck-Institut für Eisenforschung GmbH > Max-Planck-Straße 1 > D-40237 Düsseldorf > > Handelsregister B 2533 > Amtsgericht Düsseldorf > > Geschäftsführung > Prof. Dr. Gerhard Dehm > Prof. Dr. Jörg Neugebauer > Prof. Dr. Dierk Raabe > Dr. Kai de Weldige > > Ust.-Id.-Nr.: DE 11 93 58 514 > Steuernummer: 105 5891 1000 > > > Please consider that invitations and e-mails of our institute are > only valid if they end with …(a)mpie.de. > If you are not sure of the validity please contact rco(a)mpie.de <mailto:[email protected]> > > Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails > aus unserem Haus nur mit der Endung …(a)mpie.de gültig sind. > In Zweifelsfällen wenden Sie sich bitte an rco(a)mpie.de <mailto:[email protected]> > -------------------------------------------------

2 7

Question about petsc4py with cuda
by MIA 15 Jan '24

15 Jan '24

Hi, I am trying to solve a large linear equation, which needs a GPU solver as comparison. I install a CUDA-enabled PETSc and petsc4py from sources using the release tarball. According to the test results after installation, the PETSc can successfully work with cuda. All my programs are written in python, so I turn to petsc4py. But I do not find any commands that define variables on coda device or define where the direct solver is executed. I check `nvidia-smi` and find my cuda does not work at all when executing my python script: from petsc4py import PETSc import numpy as np n = 1000 nnz = 3 * np.ones(n, dtype=np.int32) nnz[0] = nnz[-1] = 2 A = PETSc.Mat() A.createAIJ([n, n], nnz=nnz) # First set the first row A.setValue(0, 0, 2) A.setValue(0, 1, -1) # Now we fill the last row A.setValue(999, 998, -1) A.setValue(999, 999, 2) # And now everything else for index in range(1, n - 1): A.setValue(index, index - 1, -1) A.setValue(index, index, 2) A.setValue(index, index + 1, -1) A.assemble() indexptr, indices, data = A.getValuesCSR() b = A.createVecLeft() b.array[:] = 1 for i in range(10): ksp = PETSc.KSP().create() ksp.setOperators(A) ksp.setType('preonly') ksp.setConvergenceHistory() ksp.getPC().setType('lu') x = A.createVecRight() ksp.solve(2*b, x) residual = A * x - 2*b if i % 10 == 0: print(f"The relative residual is: {residual.norm() / b.norm()}.”) What should I do to utlize GPU to execute the KSP task? Are there some settings to be modified? Looking forward to your early reply. Thanks a lot.

2 1

Question about petsc4py with cuda
by 1143418754 15 Jan '24

15 Jan '24

1143418754 1143418754(a)qq.com  Hi, I am trying to solve a large linear equation, which needs a GPU solver as comparison. I install a CUDA-enabled PETSc and petsc4py from sources using the release tarball. According to the test results after installation, the PETSc can successfully work with cuda.  All my programs are written in python, so I turn to petsc4py. But I do not find any commands that define variables on coda device or define where the direct solver is executed. I check `nvidia-smi` and find my cuda does not work at all when executing my python script: from petsc4py import PETSc import numpy as np n = 1000 nnz = 3 * np.ones(n, dtype=np.int32) nnz[0] = nnz[-1] = 2 A = PETSc.Mat() A.createAIJ([n, n], nnz=nnz) # First set the first row A.setValue(0, 0, 2) A.setValue(0, 1, -1) # Now we fill the last row A.setValue(999, 998, -1) A.setValue(999, 999, 2) # And now everything else for index in range(1, n - 1): A.setValue(index, index - 1, -1) A.setValue(index, index, 2) A.setValue(index, index + 1, -1)  A.assemble() indexptr, indices, data = A.getValuesCSR() b = A.createVecLeft() b.array[:] = 1 for i in range(10): ksp = PETSc.KSP().create() ksp.setOperators(A) ksp.setType('preonly') ksp.setConvergenceHistory() ksp.getPC().setType('lu') x = A.createVecRight() ksp.solve(2*b, x) residual = A * x - 2*b if i % 10 == 0: print(f"The relative residual is: {residual.norm() / b.norm()}.”) What should I do to utlize GPU to execute the KSP task? Are there some settings to be modified? Looking forward to your early reply. Thanks a lot.

1 0

Neumann Boundary Condition with DMDACreate3D
by Gourav Kumbhojkar 14 Jan '24

14 Jan '24

Hi, I am trying to implement a No-flux boundary condition for a 3D domain. I previously modeled a no flux boundary in 2D domain using DMDACreate2D and “PETSC_BOUNDARY_MIRROR” which worked well. However, the manual pages say that the Mirror Boundary is not supported for 3D. Could you please point me to the right resources to implement no flux boundary condition in 3D domains? Regards, Gourav K.

2 8

Parallel processes run significantly slower
by Steffen Wilksen | Universitaet Bremen 12 Jan '24

12 Jan '24

Hi all, I'm trying to do repeated matrix-vector-multiplication of large sparse matrices in python using petsc4py. Even the most simple method of parallelization, dividing up the calculation to run on multiple processes indenpendtly, does not seem to give a singnificant speed up for large matrices. I constructed a minimal working example, which I run using mpiexec -n N python parallel_example.py, where N is the number of processes. Instead of taking approximately the same time irrespective of the number of processes used, the calculation is much slower when starting more MPI processes. This translates to little to no speed up when splitting up a fixed number of calculations over N processes. As an example, running with N=1 takes 9s, while running with N=4 takes 34s. When running with smaller matrices, the problem is not as severe (only slower by a factor of 1.5 when setting MATSIZE=1e+5 instead of MATSIZE=1e+6). I get the same problems when just starting the script four times manually without using MPI. I attached both the script and the log file for running the script with N=4. Any help would be greatly appreciated. Calculations are done on my laptop, arch linux version 6.6.8 and PETSc version 3.20.2. Kind Regards Steffen

4 6