July 2013 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] Improving and stabilizing GPU support
by Jed Brown 22 Jul '13

22 Jul '13

Steven Dalton <sdalton1(a)gmail.com> writes: > Since CUSP provides support for wrapping generic GPU memory > into CUSP vectors it shouldn't be difficult to interface with the CUSP > preconditioners either way. Do you intend for CUSP to ever support distributed memory, or will it always be purely local?

1 0

Re: [petsc-dev] windows fixes.
by Jed Brown 21 Jul '13

21 Jul '13

Jed Brown <jedbrown(a)mcs.anl.gov> writes: > I still don't understand the actual problem. I reproduced this funny > thing where Python prints 'gcc' and a stack trace without stating an > exception, but I didn't figure out what the exception was or why the > reporting was messed up. It turns out that user dict behavior was only added in Python-2.5, so __missing__ was not invoked as expected. Pushed the fix to 'next'. Dave, can you try using 'next'?

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Steven Dalton 21 Jul '13

21 Jul '13

Hey Jed, Just subscribed. Interesting thread. It seems reasonable to provide a CUDA backend. Providing a interface around CUBLAS and/or CUSPARSE should be sufficient to provide base functionality. Since CUSP provides support for wrapping generic GPU memory into CUSP vectors it shouldn't be difficult to interface with the CUSP preconditioners either way. Thrust seems most appropriate to ease development associated with constructing nonstandard operations and/or implementing custom datatypes on the GPU. I'm still somewhat new so it's not clear to me how far Petsc will expand into those areas on the GPU. Steve On 7/20/13 1:40 PM, "Jed Brown" <jedbrown(a)mcs.anl.gov> wrote: >Steve, are you subscribed to the petsc-dev mailing list? > >Karl Rupp <rupp(a)mcs.anl.gov> writes: > >> Hi Paul, >> >> >>> * Reduce CUSP dependency: The current elementary operations are >>>>> mainly realized via CUSP. With better support via CUSPARSE and >>>>> CUBLAS, I'd add a separate 'native' CUDA backend so that we can >>>>> provide a full set of vector and sparse matrix operations out of the >>>>> default NVIDIA toolchain. We will still keep CUSP for its >>>>> preconditioners, yet we no longer depend on it. >>> Agreed. In the past, I've suggested a -vec_type cuda (not cusp). All >>>the >>> CUSP operations can be done with Thrust algorithms. Since Thrust comes >>> default with CUDA, one can have only a CUDA dependency. >> >> Yes, I opt for >> -vec_type cuda >> if everything needed is shipped with the CUDA toolkit. I even tend to >> avoid Thrust as much as possible and go with CUBLAS/CUSPARSE because we >> get faster compilation and less compiler warnings this way, but that's >> an implementation detail :-) >> >> >>>>> * Integrate last bits of txpetscgpu package. I assume Paul will >>>>> provide a helping hand here. >>> Of course. This will go much faster as much of the hard work is done. >>>Do >>> people want support for different matrix formats in the CUSP classes : >>> i.e. diagonal, ellpack, hybrid? I think the CUSP preconditioners can be >>> derived from matrices stored in non-csr format (although they're likely >>> just doing a convert under the hood). >> >> Since people keep asking for fast SpMV, we should provide these other >> formats as well (actually, they are partially provided with your update >> to the CUSPARSE bindings already). The main reason for CUSP is the SA >> preconditioner, for which SpMV performance doesn't really matter. > >Well, SpMV affects cycle time, but setup is primarily sparse >matrix-matrix. > >>>>> * Documentation: Add a chapter on GPUs to the manual, particularly on >>>>> what to expect and what not to expect. Update documentation on >>>>> webpage regarding installation. >>> I will help with the manual. >> >> Cheers :-) >> >> >>>>> * Integration of FEM quadrature from SNES ex52. The CUDA part >>>>> requiring code generation is not very elegant, while the OpenCL >>>>> approach is better suited for a library integration thanks to JIT. >>>>> However, this requires user code to be provided as a string (again >>>>> not very elegant) or loaded from file (more reasonable). How much FEM >>>>> functionality do we want to provide via PETSc? >>> Multi-GPU is a highly pressing need, IMO. Need to figure out how to >>>make >>> Block Jacobi and ASM run efficiently. >> >> The tricky part here is to balance processes vs. threads vs. GPUs. If >>we >> use more than one GPU per process, we will duplicate more and more of >> the current MPI logic over time just to move data between GPUs. >>However, >> if we just use one GPU per process, we will under-utilize the CPU >>unless >> we have a good interaction with threadcomm. >> >> Best regards, >> Karli

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Jed Brown 20 Jul '13

20 Jul '13

Karl Rupp <rupp(a)mcs.anl.gov> writes: > Jed gave me a good hint with respect to D-ILU0, which I'll also add to > PETSc. As with other GPU-accelerations using ILU, it will require a > proper matrix ordering to give good performance. Well, factorization time is not the bottleneck in ILU(0) so there's not much sense making the preconditioner weaker to get more parallelism in that part. It might be nice to have a CPU implementation just so we can compare strength. > I'm somewhat tempted to port the SA-AMG implementation in CUSP to > OpenCL as well, but this certainly won't be in 1.5.0. Rather than "porting" SACUSP to distributed memory, we are taking a different approach: providing the matrix-matrix operations for use with PCGAMG.

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Jed Brown 20 Jul '13

20 Jul '13

Steve, are you subscribed to the petsc-dev mailing list? Karl Rupp <rupp(a)mcs.anl.gov> writes: > Hi Paul, > > >>> * Reduce CUSP dependency: The current elementary operations are >>>> mainly realized via CUSP. With better support via CUSPARSE and >>>> CUBLAS, I'd add a separate 'native' CUDA backend so that we can >>>> provide a full set of vector and sparse matrix operations out of the >>>> default NVIDIA toolchain. We will still keep CUSP for its >>>> preconditioners, yet we no longer depend on it. >> Agreed. In the past, I've suggested a -vec_type cuda (not cusp). All the >> CUSP operations can be done with Thrust algorithms. Since Thrust comes >> default with CUDA, one can have only a CUDA dependency. > > Yes, I opt for > -vec_type cuda > if everything needed is shipped with the CUDA toolkit. I even tend to > avoid Thrust as much as possible and go with CUBLAS/CUSPARSE because we > get faster compilation and less compiler warnings this way, but that's > an implementation detail :-) > > >>>> * Integrate last bits of txpetscgpu package. I assume Paul will >>>> provide a helping hand here. >> Of course. This will go much faster as much of the hard work is done. Do >> people want support for different matrix formats in the CUSP classes : >> i.e. diagonal, ellpack, hybrid? I think the CUSP preconditioners can be >> derived from matrices stored in non-csr format (although they're likely >> just doing a convert under the hood). > > Since people keep asking for fast SpMV, we should provide these other > formats as well (actually, they are partially provided with your update > to the CUSPARSE bindings already). The main reason for CUSP is the SA > preconditioner, for which SpMV performance doesn't really matter. Well, SpMV affects cycle time, but setup is primarily sparse matrix-matrix. >>>> * Documentation: Add a chapter on GPUs to the manual, particularly on >>>> what to expect and what not to expect. Update documentation on >>>> webpage regarding installation. >> I will help with the manual. > > Cheers :-) > > >>>> * Integration of FEM quadrature from SNES ex52. The CUDA part >>>> requiring code generation is not very elegant, while the OpenCL >>>> approach is better suited for a library integration thanks to JIT. >>>> However, this requires user code to be provided as a string (again >>>> not very elegant) or loaded from file (more reasonable). How much FEM >>>> functionality do we want to provide via PETSc? >> Multi-GPU is a highly pressing need, IMO. Need to figure out how to make >> Block Jacobi and ASM run efficiently. > > The tricky part here is to balance processes vs. threads vs. GPUs. If we > use more than one GPU per process, we will duplicate more and more of > the current MPI logic over time just to move data between GPUs. However, > if we just use one GPU per process, we will under-utilize the CPU unless > we have a good interaction with threadcomm. > > Best regards, > Karli

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Karl Rupp 19 Jul '13

19 Jul '13

Hi, > Do you know anything about nvamg? My understanding is that it is an Nvidia > project that has not been advertised much and was used in some special apps > like Ansys Fluent. I was at a talk by an NVIDIA guy in May. It's an aggregation-based algebraic multigrid scheme, using similar components like what is implemented in CUSP. It will be a commercial product from NVIDIA (maybe free or cheap academic licenses) using MPI on a small number of nodes, primarily targeting widespread engineering software. The team consists of about 5-10 members and is headed by Robert Strzodka if I remember correctly, so they need to generate enough revenue to 'feed' them :-) > Thanks for sharing your experiences. I found the paper below earlier this > evening when I revisited the ViennaCL web site and looked at the benchmark > section. Ah, I see. :-) Even though I spend some time on improving MIC performance since these benchmarks were run, the overall trend remained unchanged. I'm still hoping for a better OpenCL SDK for MIC from Intel... Best regards, Karli

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Dave Nystrom 19 Jul '13

19 Jul '13

Karl Rupp writes: > Hi Dave, > > > That sounds very reasonable. Regarding polynomial preconditioning, were you > > thinking of least squares polynomial preconditioning or something else? > > I haven't thought about anything specific yet, just about the > infrastructure for applying any p(A). Okay. I knew that Paul had implemented a least squares polynomial preconditioner and published some results with that. So I had wondered if you might be working with that. > > > > Will there be any improvements for GPU preconditioners in ViennaCL 1.5.0? > > > > When do you expect ViennaCL 1.5.0 to be available in PETSc? > > > > > > Jed gave me a good hint with respect to D-ILU0, which I'll also add to > > > PETSc. As with other GPU-accelerations using ILU, it will require a > > > proper matrix ordering to give good performance. I'm somewhat tempted to > > > port the SA-AMG implementation in CUSP to OpenCL as well, but this > > > certainly won't be in 1.5.0. > > > > Porting SA-AMG to OpenCL also sounds attractive. I was thinking that the > > ViennaCL documentation already mentioned an algebraic preconditioner that was > > in alpha or beta status. > > The current AMG implementations all require a CPU-based setup stage and > thus limit the gain you could eventually get. In some cases where the > setup is less pronounced (e.g. lagging the preconditioner for nonlinear > or time-dependent problems) this is fine, but for stationary linear > problems with regular operators this is not very competitive. Do you know anything about nvamg? My understanding is that it is an Nvidia project that has not been advertised much and was used in some special apps like Ansys Fluent. > > I'm still trying to get my mind around the memory bandwidth issue for sparse > > linear algebra. Your report above of the Intel result adds to my confusion. > > From my understanding, the theoretical peak memory bandwidth for some systems > > of interest is as follows: > > > > Dual socket Sandy Bridge: 102 GB/s > > Nvidia Kepler K20X: 250 GB/s > > Intel Xeon Phi: 350 GB/s > > > > What I am trying to understand is what sort of memory bandwidth is achievable > > by a good implementation for the sparse linear algebra that PETSc does with > > an iterative solver like CG using Jacobi preconditioning. The plots which I > > sent links to yesterday seemed to show memory bandwidth for a dual socket > > Sandy Bridge to be well below the theoretical peak, perhaps less than 50 GB/s > > for 16 threads. For Xeon Phi, you are saying that Intel could not get more > > than 95 GB/s. But I saw a presentation last week where Nvidia was getting > > about 200 GB/s for a matrix transpose. So it makes me wonder if the > > different systems are equally good at exploiting their theoretical peak > > memory bandwidths or whether one, like the Nvidia K20X, might be better. If > > that were the case, then I might expect a good implementation of sparse > > linear algebra on a Kepler K20X to be 4-5 times faster than a good > > implementation on a dual socket Sandy Bridge node rather than a 2.5x > > difference. > > Intel's marketing machinery was tricking you: The 350 GB/sec are the > peak bandwidth from the ring bus connecting the MIC cores to GDDRAM. > However, the internal ring bus operates at only 220 GB/sec (see for > example the following paper [1]). With some prefetching tricks and Intel > pragma/compiler magic one obtains about 160 GB/sec for the STREAM > benchmark, which is 75% of peak. The Intel OpenCL SDK adds another loss > here, resulting in only 95 GB/sec. This was why I got in contact with > Intel in order to find out whether this is a weakness of the SDK or > whether I missed something. Turned out to be the former... > > As you know, for dual Socket systems one only gets good bandwidth if the > placement in memory is done in order to adhere to NUMA. On such a dual > socket system I recently managed to get 75 GB/sec with OpenCL, which is > again 75% of peak performance. Unfortunately OpenCL does not consider > NUMA, so this is not very stable, so you may only get half of it if all > data happens to reside on the same memory link. > > On GPUs including the K20X one also obtains about 75% of peak: On a > Radeon 7970 I got 220 out of 288 theoretical peak, other people even > reported up to 250 GB/sec for a GTX Titan (288 GB/sec theoretical peak), > and I also got 131 GB/sec out of 159 GB/sec peak for a rather dated GTX 285. > > Overall, the rule of thumb seems to be 75% of peak if everything is done > correctly and if one finds the right baseline (the Xeon Phi is a beast > in this regard). These are numbers for sequential reads, hence no cache > effects or other mechanisms such as paging cause other spurious effects. > When it comes to actual optimizations for sparse linear algebra, CPUs > and GPUs ask for slightly different sets of optimizations because cache > lines and memory controllers differ... Thanks for sharing your experiences. I found the paper below earlier this evening when I revisited the ViennaCL web site and looked at the benchmark section. > Best regards, > > Karli > > [1] http://arxiv.org/abs/1302.1078

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Karl Rupp 19 Jul '13

19 Jul '13

Hi Dave, > That sounds very reasonable. Regarding polynomial preconditioning, were you > thinking of least squares polynomial preconditioning or something else? I haven't thought about anything specific yet, just about the infrastructure for applying any p(A). > > > Will there be any improvements for GPU preconditioners in ViennaCL 1.5.0? > > > When do you expect ViennaCL 1.5.0 to be available in PETSc? > > > > Jed gave me a good hint with respect to D-ILU0, which I'll also add to > > PETSc. As with other GPU-accelerations using ILU, it will require a > > proper matrix ordering to give good performance. I'm somewhat tempted to > > port the SA-AMG implementation in CUSP to OpenCL as well, but this > > certainly won't be in 1.5.0. > > Porting SA-AMG to OpenCL also sounds attractive. I was thinking that the > ViennaCL documentation already mentioned an algebraic preconditioner that was > in alpha or beta status. The current AMG implementations all require a CPU-based setup stage and thus limit the gain you could eventually get. In some cases where the setup is less pronounced (e.g. lagging the preconditioner for nonlinear or time-dependent problems) this is fine, but for stationary linear problems with regular operators this is not very competitive. > I'm still trying to get my mind around the memory bandwidth issue for sparse > linear algebra. Your report above of the Intel result adds to my confusion. > From my understanding, the theoretical peak memory bandwidth for some systems > of interest is as follows: > > Dual socket Sandy Bridge: 102 GB/s > Nvidia Kepler K20X: 250 GB/s > Intel Xeon Phi: 350 GB/s > > What I am trying to understand is what sort of memory bandwidth is achievable > by a good implementation for the sparse linear algebra that PETSc does with > an iterative solver like CG using Jacobi preconditioning. The plots which I > sent links to yesterday seemed to show memory bandwidth for a dual socket > Sandy Bridge to be well below the theoretical peak, perhaps less than 50 GB/s > for 16 threads. For Xeon Phi, you are saying that Intel could not get more > than 95 GB/s. But I saw a presentation last week where Nvidia was getting > about 200 GB/s for a matrix transpose. So it makes me wonder if the > different systems are equally good at exploiting their theoretical peak > memory bandwidths or whether one, like the Nvidia K20X, might be better. If > that were the case, then I might expect a good implementation of sparse > linear algebra on a Kepler K20X to be 4-5 times faster than a good > implementation on a dual socket Sandy Bridge node rather than a 2.5x > difference. Intel's marketing machinery was tricking you: The 350 GB/sec are the peak bandwidth from the ring bus connecting the MIC cores to GDDRAM. However, the internal ring bus operates at only 220 GB/sec (see for example the following paper [1]). With some prefetching tricks and Intel pragma/compiler magic one obtains about 160 GB/sec for the STREAM benchmark, which is 75% of peak. The Intel OpenCL SDK adds another loss here, resulting in only 95 GB/sec. This was why I got in contact with Intel in order to find out whether this is a weakness of the SDK or whether I missed something. Turned out to be the former... As you know, for dual Socket systems one only gets good bandwidth if the placement in memory is done in order to adhere to NUMA. On such a dual socket system I recently managed to get 75 GB/sec with OpenCL, which is again 75% of peak performance. Unfortunately OpenCL does not consider NUMA, so this is not very stable, so you may only get half of it if all data happens to reside on the same memory link. On GPUs including the K20X one also obtains about 75% of peak: On a Radeon 7970 I got 220 out of 288 theoretical peak, other people even reported up to 250 GB/sec for a GTX Titan (288 GB/sec theoretical peak), and I also got 131 GB/sec out of 159 GB/sec peak for a rather dated GTX 285. Overall, the rule of thumb seems to be 75% of peak if everything is done correctly and if one finds the right baseline (the Xeon Phi is a beast in this regard). These are numbers for sequential reads, hence no cache effects or other mechanisms such as paging cause other spurious effects. When it comes to actual optimizations for sparse linear algebra, CPUs and GPUs ask for slightly different sets of optimizations because cache lines and memory controllers differ... Best regards, Karli [1] http://arxiv.org/abs/1302.1078

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Dave Nystrom 19 Jul '13

19 Jul '13

Hi Karli, Karl Rupp writes: > Hi Dave, > > > Your list sounds great to me. Glad that you and Paul are working on > > this together. > > > > My main interests are in better preconditioner support and better > > multi-GPU/MPI scalability. > > This is follow-up work then. There are a couple of 'simple' > preconditioners (polynomial preconditioning, maybe some point-block > Jacobi) which can also be useful as smoothers and which we can add in the > near future. We should just get the 'infrastructure' work done first so > that we don't have to unnecessarily adjust too much code later on. That sounds very reasonable. Regarding polynomial preconditioning, were you thinking of least squares polynomial preconditioning or something else? > > Is there any progress on Steve Dalton's work on the cusp algebraic multigrid > > preconditioner with PETSc? I believe Jed said in a previous email that Steve > > was going to be working on adding MPI support for that as well as other > > enhancements. > > Yes, Steve is working on this right here at our division. Jed can give a > more detailed answer on this. Glad to hear that that project is going forward. > > Will there be any improvements for GPU preconditioners in ViennaCL 1.5.0? > > When do you expect ViennaCL 1.5.0 to be available in PETSc? > > Jed gave me a good hint with respect to D-ILU0, which I'll also add to > PETSc. As with other GPU-accelerations using ILU, it will require a > proper matrix ordering to give good performance. I'm somewhat tempted to > port the SA-AMG implementation in CUSP to OpenCL as well, but this > certainly won't be in 1.5.0. Porting SA-AMG to OpenCL also sounds attractive. I was thinking that the ViennaCL documentation already mentioned an algebraic preconditioner that was in alpha or beta status. > > I'm also interested in trying the PETSc ViennaCL support on the Xeon Phi. > > Do you have a schedule for when that might be ready for friendly testing? > > With OpenCL you can already test this now. Just install the Intel OpenCL > SDK on your Xeon Phi machine, configure with --download-viennacl, > --with-opencl-include=..., --with-opencl-lib=..., and pass the > -viennacl_device_accelerator > flag in addition to -vec_type viennacl -mat_type aijviennacl when executing. > > Unfortunately the application memory bandwidth we get on the Xeon Phi is > too limited to be useful for off-loaded execution as it is the case with > OpenCL: Even the folks at Intel couldn't obtain more than ~95 GB/sec > even when filling up the whole MIC with just two vectors for > benchmarking a simple copy operation. I'm still trying to get my mind around the memory bandwidth issue for sparse linear algebra. Your report above of the Intel result adds to my confusion. >From my understanding, the theoretical peak memory bandwidth for some systems of interest is as follows: Dual socket Sandy Bridge: 102 GB/s Nvidia Kepler K20X: 250 GB/s Intel Xeon Phi: 350 GB/s What I am trying to understand is what sort of memory bandwidth is achievable by a good implementation for the sparse linear algebra that PETSc does with an iterative solver like CG using Jacobi preconditioning. The plots which I sent links to yesterday seemed to show memory bandwidth for a dual socket Sandy Bridge to be well below the theoretical peak, perhaps less than 50 GB/s for 16 threads. For Xeon Phi, you are saying that Intel could not get more than 95 GB/s. But I saw a presentation last week where Nvidia was getting about 200 GB/s for a matrix transpose. So it makes me wonder if the different systems are equally good at exploiting their theoretical peak memory bandwidths or whether one, like the Nvidia K20X, might be better. If that were the case, then I might expect a good implementation of sparse linear algebra on a Kepler K20X to be 4-5 times faster than a good implementation on a dual socket Sandy Bridge node rather than a 2.5x difference. > Thus, I don't think our efforts are currently well spent on trying a fully > native execution of PETSc on the MIC, because the trend is going more > towards a tighter CPU/accelerator integration on the same die rather than > piggy-backing via PCI-Express. Anyway, I'll let you know if there are any > updates on this front. Thanks. Please do. Best regards, Dave > Best regards, > > Karli

1 0

Re: [petsc-dev] Improving and stabilizing GPU support
by Karl Rupp 19 Jul '13

19 Jul '13

Hi Paul, >>> * Reduce CUSP dependency: The current elementary operations are >>> mainly realized via CUSP. With better support via CUSPARSE and >>> CUBLAS, I'd add a separate 'native' CUDA backend so that we can >>> provide a full set of vector and sparse matrix operations out of the >>> default NVIDIA toolchain. We will still keep CUSP for its >>> preconditioners, yet we no longer depend on it. > Agreed. In the past, I've suggested a -vec_type cuda (not cusp). All the > CUSP operations can be done with Thrust algorithms. Since Thrust comes > default with CUDA, one can have only a CUDA dependency. Yes, I opt for -vec_type cuda if everything needed is shipped with the CUDA toolkit. I even tend to avoid Thrust as much as possible and go with CUBLAS/CUSPARSE because we get faster compilation and less compiler warnings this way, but that's an implementation detail :-) >>> * Integrate last bits of txpetscgpu package. I assume Paul will >>> provide a helping hand here. > Of course. This will go much faster as much of the hard work is done. Do > people want support for different matrix formats in the CUSP classes : > i.e. diagonal, ellpack, hybrid? I think the CUSP preconditioners can be > derived from matrices stored in non-csr format (although they're likely > just doing a convert under the hood). Since people keep asking for fast SpMV, we should provide these other formats as well (actually, they are partially provided with your update to the CUSPARSE bindings already). The main reason for CUSP is the SA preconditioner, for which SpMV performance doesn't really matter. >>> * Documentation: Add a chapter on GPUs to the manual, particularly on >>> what to expect and what not to expect. Update documentation on >>> webpage regarding installation. > I will help with the manual. Cheers :-) >>> * Integration of FEM quadrature from SNES ex52. The CUDA part >>> requiring code generation is not very elegant, while the OpenCL >>> approach is better suited for a library integration thanks to JIT. >>> However, this requires user code to be provided as a string (again >>> not very elegant) or loaded from file (more reasonable). How much FEM >>> functionality do we want to provide via PETSc? > Multi-GPU is a highly pressing need, IMO. Need to figure out how to make > Block Jacobi and ASM run efficiently. The tricky part here is to balance processes vs. threads vs. GPUs. If we use more than one GPU per process, we will duplicate more and more of the current MPI logic over time just to move data between GPUs. However, if we just use one GPU per process, we will under-utilize the CPU unless we have a good interaction with threadcomm. Best regards, Karli

1 0