<div dir="ltr">I do one complete solve to get everything setup, to be safe.<div><br></div><div>src/ts/tutorials/ex13.c does this and runs multiple solves, if you like but one solve is probably fine.</div><div>This was designed as a benchmark and is nice because it can do any order FE solve of Poisson (uses DM/PetscFE, slow).</div><div>src/ksp/ksp/tutorials/ex56.c is old school, hardwired for elasticity but is simpler and the setup is faster if you are doing large problems per MPI process.<br></div><div><br></div><div>Mark</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Feb 6, 2023 at 2:06 PM Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div><br></div> It should not crash, take a look at the test cases at the bottom of the file. You are likely correct if the code, unfortunately, does use DMCreateMatrix() it will not work out of the box for geometric multigrid. So it might be the wrong example for you.<div><br></div><div>  I don't know what you mean about clever. If you simply set the solution to zero at the beginning of the loop then it will just do the same solve multiple times. The setup should not do much of anything after the first solver.  Thought usually solves are big enough that one need not run solves multiple times to get a good understanding of their performance.</div><div><br></div><div><br></div><div>  <br><div><br></div><div><br><div><br><blockquote type="cite"><div>On Feb 6, 2023, at 12:44 PM, Paul Grosse-Bley <<a href="mailto:paul.grosse-bley@ziti.uni-heidelberg.de" target="_blank">paul.grosse-bley@ziti.uni-heidelberg.de</a>> wrote:</div><br><div>Hi Barry,<br><br>src/ksp/ksp/tutorials/bench_kspsolve.c is certainly the better starting point, thank you! Sadly I get a segfault when executing that example with PCMG and more than one level, i.e. with the minimal args:<br><br>$ mpiexec -c 1 ./bench_kspsolve -split_ksp -pc_type mg -pc_mg_levels 2<br>===========================================<br>Test: KSP performance - Poisson<br>    Input matrix: 27-pt finite difference stencil<br>    -n 100<br>    DoFs = 1000000<br>    Number of nonzeros = 26463592<br><br>Step1  - creating Vecs and Mat...<br>Step2a - running PCSetUp()...<br>[0]PETSC ERROR: ------------------------------------------------------------------------<br>[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range<br>[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger<br>[0]PETSC ERROR: or see <a href="https://petsc.org/release/faq/#valgrind" target="_blank">https://petsc.org/release/faq/#valgrind</a> and <a href="https://petsc.org/release/faq/" target="_blank">https://petsc.org/release/faq/</a><br>[0]PETSC ERROR: or try <a href="https://docs.nvidia.com/cuda/cuda-memcheck/index.html" target="_blank">https://docs.nvidia.com/cuda/cuda-memcheck/index.html</a> on NVIDIA CUDA systems to find memory corruption errors<br>[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run<br>[0]PETSC ERROR: to get more information on the crash.<br>[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.<br>--------------------------------------------------------------------------<br>MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD<br>with errorcode 59.<br><br>NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.<br>You may or may not see output from other processes, depending on<br>exactly when Open MPI kills them.<br>--------------------------------------------------------------------------<br><br>As the matrix is not created using DMDACreate3d I expected it to fail due to the missing geometric information, but I expected it to fail more gracefully than with a segfault.<br>I will try to combine bench_kspsolve.c with ex45.c to get easy MG preconditioning, especially since I am interested in the 7pt stencil for now.<br><br>Concerning my benchmarking loop from before: Is it generally discouraged to do this for KSPSolve due to PETSc cleverly/lazily skipping some of the work when doing the same solve multiple times or are the solves not iterated in bench_kspsolve.c (while the MatMuls are with -matmult) just to keep the runtime short?<br><br>Thanks,<br>Paul<br><br>On Monday, February 06, 2023 17:04 CET, Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br> <blockquote type="cite"> </blockquote><div> </div>  Paul,<div> </div><div>   I think src/ksp/ksp/tutorials/benchmark_ksp.c is the code intended to be used for simple benchmarking. </div><div> </div><div>   You can use VecCudaGetArray() to access the GPU memory of the vector and then call a CUDA kernel to compute the right hand side vector directly on the GPU.</div><div> </div><div>  Barry</div><div> <div> <blockquote type="cite"><div>On Feb 6, 2023, at 10:57 AM, Paul Grosse-Bley <<a href="mailto:paul.grosse-bley@ziti.uni-heidelberg.de" target="_blank">paul.grosse-bley@ziti.uni-heidelberg.de</a>> wrote:</div> <div>Hi,<br><br>I want to compare different implementations of multigrid solvers for Nvidia GPUs using the poisson problem (starting from ksp tutorial example ex45.c).<br>Therefore I am trying to get runtime results comparable to <a href="https://bitbucket.org/nsakharnykh/hpgmg-cuda/src/master/" target="_blank">hpgmg-cuda</a> (finite-volume), i.e. using multiple warmup and measurement solves and avoiding measuring setup time.<br>For now I am using -log_view with added stages:<br><br>PetscLogStageRegister("Solve Bench", &solve_bench_stage);<br>  for (int i = 0; i < BENCH_SOLVES; i++) {<br>    PetscCall(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, NULL)); // reset x<br>    PetscCall(KSPSetUp(ksp)); // try to avoid setup overhead during solve<br>    PetscCall(PetscDeviceContextSynchronize(dctx)); // make sure that everything is done<br><br>    PetscLogStagePush(solve_bench_stage);<br>    PetscCall(KSPSolve(ksp, NULL, NULL));<br>    PetscLogStagePop();<br>  }<br><br>This snippet is preceded by a similar loop for warmup.<br><br>When profiling this using Nsight Systems, I see that the very first solve is much slower which mostly correspods to H2D (host to device) copies and e.g. cuBLAS setup (maybe also paging overheads as mentioned in the <a href="https://petsc.org/release/docs/manual/profiling/#accurate-profiling-and-paging-overheads" target="_blank">docs</a>, but probably insignificant in this case). The following solves have some overhead at the start from a H2D copy of a vector (the RHS I guess, as the copy is preceeded by a matrix-vector product) in the first MatResidual call (callchain: KSPSolve->MatResidual->VecAYPX->VecCUDACopyTo->cudaMemcpyAsync). My interpretation of the profiling results (i.e. cuBLAS calls) is that that vector is overwritten with the residual in Daxpy and therefore has to be copied again for the next iteration.<br><br>Is there an elegant way of avoiding that H2D copy? I have seen some examples on constructing matrices directly on the GPU, but nothing about vectors. Any further tips for benchmarking (vs profiling) PETSc solvers? At the moment I am using jacobi as smoother, but I would like to have a CUDA implementation of SOR instead. Is there a good way of achieving that, e.g. using PCHYPREs boomeramg with a single level and "SOR/Jacobi"-smoother  as smoother in PCMG? Or is the overhead from constantly switching between PETSc and hypre too big?<br><br>Thanks,<br>Paul</div></blockquote></div></div>
</div></blockquote></div><br></div></div></div></blockquote></div>