February 2016 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Jose E. Roman 28 Feb '16

28 Feb '16

> El 28 feb 2016, a las 10:45, Karl Rupp <rupp(a)iue.tuwien.ac.at> escribió: > > Hi, > >> I like the idea of having separate VECCUDA and VECVIENNACL, because it is possible to implement VECCUDA without dependence on a C++ compiler (only the CUDA compiler). > > I don't understand this part. NVCC also requires a C++ host compiler and is fairly picky about the supported compilers. You are right. I was thinking of the case when one has a pure C code and wants to use a --with-language=C PETSc configuration. > > >> If you want, we can prepare a rough initial implementation of VECCUDA in the next days, and we can later discuss what to keep/discard. > > Any contributions are welcome :-) > > >> Karl: regarding the time constraints, our idea is to present something at a conference this summer, and deadlines are approaching. > > Ok, this is on fairly short notice considering the changes required. I recommend to start with copying the CUSP sources and migrate it over to VECCUDA by replacing any use of cusp::array1d to a raw CUDA handle. Operations from CUSP should be replaced by CUBLAS calls. Ok. Will start work on this. Jose > > Best regards, > Karli >

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Karl Rupp 28 Feb '16

28 Feb '16

Hi, > I like the idea of having separate VECCUDA and VECVIENNACL, because it is possible to implement VECCUDA without dependence on a C++ compiler (only the CUDA compiler). I don't understand this part. NVCC also requires a C++ host compiler and is fairly picky about the supported compilers. > If you want, we can prepare a rough initial implementation of VECCUDA in the next days, and we can later discuss what to keep/discard. Any contributions are welcome :-) > Karl: regarding the time constraints, our idea is to present something at a conference this summer, and deadlines are approaching. Ok, this is on fairly short notice considering the changes required. I recommend to start with copying the CUSP sources and migrate it over to VECCUDA by replacing any use of cusp::array1d to a raw CUDA handle. Operations from CUSP should be replaced by CUBLAS calls. Best regards, Karli >>>> >>>> If there is interest we can help in adding this stuff. >>> >>> What are your time constraints? >>> >>> Best regards, >>> Karli >>> >>> >> >> -- >> Dominic Meiser >> Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303 >

1 0

Re: [petsc-dev] TSComputeRHSJacobianConstant issue
by Barry Smith 27 Feb '16

27 Feb '16

Jed and Emil, Yes, this is a terrible mess. The TSComputeRHSJacobianConstant() will only work for methods that do not change entries in the Jacobian matrix. And for some reason the runex3_4 and _5 which do not work are turned off in the makefile ! It looks like Rosw does work because it does not change the matrix? And yet it is still turned off? I don't know what a good fix is. If we want to keep this capability then the TSComputeRHSJacobianConstant() interface has to keep a copy of the original user provided matrix and copy it into the matrices used by the method each time. This may lead to ugly special case code. Something has to be done about (simplest thing is just tossing the entire concept) because now there is a very broken API in TS that beginning users are likely to try and get upset with. Barry > On Feb 26, 2016, at 11:04 AM, Dominic Meiser <dmeiser(a)txcorp.com> wrote: > > Hi all, > > I'm running into problems with TS with a time independent, linear > rhs. Following the user's manual, I use > > ierr=TSSetProblemType(ts, TS_LINEAR);CHKERRQ(ierr); > ierr=TSSetRHSFunction(ts, 0, TSComputeRHSFunctionLinear, 0);CHKERRQ(ierr); > ierr=TSSetRHSJacobian(ts, mat, mat, TSComputeRHSJacobianConstant, 0);CHKERRQ(ierr); > > But running the code I get: > > [dmeiser@ivyamd src]$ ./sd -ts_monitor -snes_monitor -ksp_monitor_true_residual -snes_converged_reason -ts_type beuler -ts_max_steps 20 TS dt 0.01 time 0. > 0 SNES Function norm 1.433485640883e+00 > 0 KSP preconditioned resid norm 1.433485640883e-02 true resid norm 1.433485640883e+00 ||r(i)||/||b|| 1.000000000000e+00 > 1 KSP preconditioned resid norm 6.528150237447e-18 true resid norm 6.122666587510e-16 ||r(i)||/||b|| 4.271173992185e-16 > 1 SNES Function norm 1.419007435910e+02 > Nonlinear solve converged due to CONVERGED_ITS iterations 1 > 1 TS dt 0.01 time 0.01 > 0 SNES Function norm 1.433342292319e+02 > Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0 > [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- > [0]PETSC ERROR: > [0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase -ts_max_snes_failures or make negative to attempt recovery > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. > [0]PETSC ERROR: Petsc Development GIT revision: pre-tsfc-330-gda8116b GIT Date: 2016-02-23 15:17:17 -0600 > [0]PETSC ERROR: ./sd on a cplx named ivyamd.txcorp.com by dmeiser Fri Feb 26 11:46:29 2016 > [0]PETSC ERROR: Configure options --with-scalar-type=complex --with-x=0 --with-debugging=1 --with-ssl=0 --with-fortran-kernels=0 --with-pthread=0 --download-fftw=1 --with-mpi=1 --download-mpich=1 --with-mpiuni-fortran-binding=0 --with-fortran-interfaces > [0]PETSC ERROR: #1 TSStep() line 3346 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c > [0]PETSC ERROR: #2 TSSolve() line 3544 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c > [0]PETSC ERROR: #3 main() line 63 in /scr_ivyamd/dmeiser/SD/src/sd.c > [0]PETSC ERROR: PETSc Option Table entries: > [0]PETSC ERROR: -ksp_monitor_true_residual > [0]PETSC ERROR: -snes_converged_reason > [0]PETSC ERROR: -snes_monitor > [0]PETSC ERROR: -ts_max_steps 2 > [0]PETSC ERROR: -ts_monitor > [0]PETSC ERROR: -ts_type beuler > [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint(a)mcs.anl.gov---------- > application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 > [unset]: aborting job: > application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 > > > It seems that the linear system is being solved just fine. But > apparently it's the wrong linear problem because the SNES > function norm doesn't decrease. > > I think the problem might be related to the following behaviour > of `scr/ts/examples/tutorials/ex3.c`: > > > [dmeiser@ivyamd tutorials]$ ./ex3 -nox -ts_type beuler -ts_max_steps 2 -ksp_converged_reason -snes_monitor > Solving a linear TS problem on 1 processor > Timestep 0: step size = 0.000143637, time = 0., 2-norm error = 1.01507e-15, max norm error = 3.10862e-15 > 0 SNES Function norm 2.018485956680e+03 > Linear solve converged due to CONVERGED_RTOL iterations 1 > 1 SNES Function norm 1.207665524289e+05 > Timestep 1: step size = 0.000143637, time = 0.000143637, 2-norm error = 0.000599907, max norm error = 0.000863771 > 0 SNES Function norm 1.195422385486e+05 > Linear solve did not converge due to DIVERGED_PCSETUP_FAILED iterations 0 > PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT > [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- > [0]PETSC ERROR: > [0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase -ts_max_snes_failures or make negative to attempt recovery > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. > [0]PETSC ERROR: Petsc Development GIT revision: pre-tsfc-330-gda8116b GIT Date: 2016-02-23 15:17:17 -0600 > [0]PETSC ERROR: ./ex3 on a cplx named ivyamd.txcorp.com by dmeiser Fri Feb 26 11:36:28 2016 > [0]PETSC ERROR: Configure options --with-scalar-type=complex --with-x=0 --with-debugging=1 --with-ssl=0 --with-fortran-kernels=0 --with-pthread=0 --download-fftw=1 --with-mpi=1 --download-mpich=1 --with-mpiuni-fortran-binding=0 --with-fortran-interfaces > [0]PETSC ERROR: #1 TSStep() line 3346 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c > [0]PETSC ERROR: #2 TSSolve() line 3544 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c > [0]PETSC ERROR: #3 main() line 216 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/examples/tutorials/ex3.c > [0]PETSC ERROR: PETSc Option Table entries: > [0]PETSC ERROR: -ksp_converged_reason > [0]PETSC ERROR: -nox > [0]PETSC ERROR: -snes_monitor > [0]PETSC ERROR: -ts_max_steps 2 > [0]PETSC ERROR: -ts_type beuler > [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint(a)mcs.anl.gov---------- > application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 > [unset]: aborting job: > application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 > > > I'm aware of Jed's commit e1244c6 which I thought was supposed to > enable constant linear rhs without IFunction and IJacobian. Note > that e1244c6 also added tests that apparently produced very > different output (see e.g. > src/ts/examples/tutorials/output/ex3_4.out). > > Is constant rhs without IFunction and IJacobian no longer > supported or did it just breat at some point? Do I need to call > TSRHSJacobianSetReuse in addition to the lines of code above? Is > it better to treat such systems using TSComputeIFunctionLinear > and TSComputeIJacobianConstant? > > Thanks, > Dominic > > -- > Dominic Meiser > Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0

PETSc release soon, request for input on needed fixes or enhancements
by Barry Smith 27 Feb '16

27 Feb '16

PETSc Users, We are planning the PETSc release 3.7 shortly. If you know of any bugs that need to be fixed or enhancements added before the release please let us know. You can think of the master branch of the PETSc repository obtainable with git clone https://bitbucket.org/petsc/petsc petsc as a release candidate for 3.7. Changes for the release are listed at http://www.mcs.anl.gov/petsc/documentation/changes/dev.html Thanks Barry

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Jose E. Roman 26 Feb '16

26 Feb '16

> El 26 feb 2016, a las 18:31, Dominic Meiser <dmeiser(a)txcorp.com> escribió: > > On Fri, Feb 26, 2016 at 02:49:39PM +0100, Karl Rupp wrote: >> >>>> The alternative would be to use raw cuda pointers instead of cusp >>>> arrays for GPU memory in VecCUSP. That would be a fairly >>>> significant undertaking (certainly more than the 2-3 weeks Karli >>>> is estimating for getting the ViennaCL cuda backend in). >>> >>> Do you mean creating a new class VECCUDA in addition to VECCUSP and VECVIENNACL? This could be a solution for us. It would mean maybe refactoring MATAIJCUSPARSE to work with these new Vecs? >> >> I prefer to replace VECCUSP with e.g. VECCUDA (and eventually also >> rename VECCUSPARSE to VECCUDA to have a unified naming for all the >> things provided natively with the CUDA SDK) in order to reduce >> external dependencies. CUSP will provide matrices, preconditioners, >> etc. as before, but is only optional and thus less likely to cause >> installation troubles. Supporting VECCUSP and VECCUDA next to each >> other is going to be too much code duplication without any benefit. > > That makes sense. At the Vec level we should be using a low > level construct (i.e. cuda raw pointers) because clients can > always provide raw pointers and they know how to consume them > (e.g. if they want to use cusp vectors on their end). > >> >> Even if we do provide VECCUDA, I still dislike the fact that we >> would have to maintain essentially the same code twice: One for >> CUDA, one for OpenCL/ViennaCL. With the ViennaCL bindings providing >> OpenMP and CUDA support soon, this also duplicates functionality. A >> possible 'fix' is to just use ViennaCL for CUDA+OpenCL+OpenMP and >> thus only maintain a single PETSc plugin for all three. However, I'm >> certainly too biased to be taken seriously here. > > I agree with this in principle. Perhaps it's time to consolidate > the cuda/cusp/cusparse/opencl efforts. Note however that > MATAIJCUSPARSE provides capabilities that won't be available > right away with ViennaCL (e.g. multi-GPU block Jacobi and ASM > preconditioners). I like the idea of having separate VECCUDA and VECVIENNACL, because it is possible to implement VECCUDA without dependence on a C++ compiler (only the CUDA compiler). If you want, we can prepare a rough initial implementation of VECCUDA in the next days, and we can later discuss what to keep/discard. Karl: regarding the time constraints, our idea is to present something at a conference this summer, and deadlines are approaching. > > Cheers, > Dominic > > >> >>> >>> If there is interest we can help in adding this stuff. >> >> What are your time constraints? >> >> Best regards, >> Karli >> >> > > -- > Dominic Meiser > Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Dominic Meiser 26 Feb '16

26 Feb '16

On Fri, Feb 26, 2016 at 02:49:39PM +0100, Karl Rupp wrote: > > >>The alternative would be to use raw cuda pointers instead of cusp > >>arrays for GPU memory in VecCUSP. That would be a fairly > >>significant undertaking (certainly more than the 2-3 weeks Karli > >>is estimating for getting the ViennaCL cuda backend in). > > > >Do you mean creating a new class VECCUDA in addition to VECCUSP and VECVIENNACL? This could be a solution for us. It would mean maybe refactoring MATAIJCUSPARSE to work with these new Vecs? > > I prefer to replace VECCUSP with e.g. VECCUDA (and eventually also > rename VECCUSPARSE to VECCUDA to have a unified naming for all the > things provided natively with the CUDA SDK) in order to reduce > external dependencies. CUSP will provide matrices, preconditioners, > etc. as before, but is only optional and thus less likely to cause > installation troubles. Supporting VECCUSP and VECCUDA next to each > other is going to be too much code duplication without any benefit. That makes sense. At the Vec level we should be using a low level construct (i.e. cuda raw pointers) because clients can always provide raw pointers and they know how to consume them (e.g. if they want to use cusp vectors on their end). > > Even if we do provide VECCUDA, I still dislike the fact that we > would have to maintain essentially the same code twice: One for > CUDA, one for OpenCL/ViennaCL. With the ViennaCL bindings providing > OpenMP and CUDA support soon, this also duplicates functionality. A > possible 'fix' is to just use ViennaCL for CUDA+OpenCL+OpenMP and > thus only maintain a single PETSc plugin for all three. However, I'm > certainly too biased to be taken seriously here. I agree with this in principle. Perhaps it's time to consolidate the cuda/cusp/cusparse/opencl efforts. Note however that MATAIJCUSPARSE provides capabilities that won't be available right away with ViennaCL (e.g. multi-GPU block Jacobi and ASM preconditioners). Cheers, Dominic > > > > > If there is interest we can help in adding this stuff. > > What are your time constraints? > > Best regards, > Karli > > -- Dominic Meiser Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0

TSComputeRHSJacobianConstant issue
by Dominic Meiser 26 Feb '16

26 Feb '16

Hi all, I'm running into problems with TS with a time independent, linear rhs. Following the user's manual, I use ierr=TSSetProblemType(ts, TS_LINEAR);CHKERRQ(ierr); ierr=TSSetRHSFunction(ts, 0, TSComputeRHSFunctionLinear, 0);CHKERRQ(ierr); ierr=TSSetRHSJacobian(ts, mat, mat, TSComputeRHSJacobianConstant, 0);CHKERRQ(ierr); But running the code I get: [dmeiser@ivyamd src]$ ./sd -ts_monitor -snes_monitor -ksp_monitor_true_residual -snes_converged_reason -ts_type beuler -ts_max_steps 20 TS dt 0.01 time 0. 0 SNES Function norm 1.433485640883e+00 0 KSP preconditioned resid norm 1.433485640883e-02 true resid norm 1.433485640883e+00 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 6.528150237447e-18 true resid norm 6.122666587510e-16 ||r(i)||/||b|| 4.271173992185e-16 1 SNES Function norm 1.419007435910e+02 Nonlinear solve converged due to CONVERGED_ITS iterations 1 1 TS dt 0.01 time 0.01 0 SNES Function norm 1.433342292319e+02 Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0 [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: [0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase -ts_max_snes_failures or make negative to attempt recovery [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Development GIT revision: pre-tsfc-330-gda8116b GIT Date: 2016-02-23 15:17:17 -0600 [0]PETSC ERROR: ./sd on a cplx named ivyamd.txcorp.com by dmeiser Fri Feb 26 11:46:29 2016 [0]PETSC ERROR: Configure options --with-scalar-type=complex --with-x=0 --with-debugging=1 --with-ssl=0 --with-fortran-kernels=0 --with-pthread=0 --download-fftw=1 --with-mpi=1 --download-mpich=1 --with-mpiuni-fortran-binding=0 --with-fortran-interfaces [0]PETSC ERROR: #1 TSStep() line 3346 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #2 TSSolve() line 3544 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #3 main() line 63 in /scr_ivyamd/dmeiser/SD/src/sd.c [0]PETSC ERROR: PETSc Option Table entries: [0]PETSC ERROR: -ksp_monitor_true_residual [0]PETSC ERROR: -snes_converged_reason [0]PETSC ERROR: -snes_monitor [0]PETSC ERROR: -ts_max_steps 2 [0]PETSC ERROR: -ts_monitor [0]PETSC ERROR: -ts_type beuler [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint(a)mcs.anl.gov---------- application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 [unset]: aborting job: application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 It seems that the linear system is being solved just fine. But apparently it's the wrong linear problem because the SNES function norm doesn't decrease. I think the problem might be related to the following behaviour of `scr/ts/examples/tutorials/ex3.c`: [dmeiser@ivyamd tutorials]$ ./ex3 -nox -ts_type beuler -ts_max_steps 2 -ksp_converged_reason -snes_monitor Solving a linear TS problem on 1 processor Timestep 0: step size = 0.000143637, time = 0., 2-norm error = 1.01507e-15, max norm error = 3.10862e-15 0 SNES Function norm 2.018485956680e+03 Linear solve converged due to CONVERGED_RTOL iterations 1 1 SNES Function norm 1.207665524289e+05 Timestep 1: step size = 0.000143637, time = 0.000143637, 2-norm error = 0.000599907, max norm error = 0.000863771 0 SNES Function norm 1.195422385486e+05 Linear solve did not converge due to DIVERGED_PCSETUP_FAILED iterations 0 PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: [0]PETSC ERROR: TSStep has failed due to DIVERGED_NONLINEAR_SOLVE, increase -ts_max_snes_failures or make negative to attempt recovery [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Development GIT revision: pre-tsfc-330-gda8116b GIT Date: 2016-02-23 15:17:17 -0600 [0]PETSC ERROR: ./ex3 on a cplx named ivyamd.txcorp.com by dmeiser Fri Feb 26 11:36:28 2016 [0]PETSC ERROR: Configure options --with-scalar-type=complex --with-x=0 --with-debugging=1 --with-ssl=0 --with-fortran-kernels=0 --with-pthread=0 --download-fftw=1 --with-mpi=1 --download-mpich=1 --with-mpiuni-fortran-binding=0 --with-fortran-interfaces [0]PETSC ERROR: #1 TSStep() line 3346 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #2 TSSolve() line 3544 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #3 main() line 216 in /scr_ivyamd/dmeiser/SD/petsc/src/ts/examples/tutorials/ex3.c [0]PETSC ERROR: PETSc Option Table entries: [0]PETSC ERROR: -ksp_converged_reason [0]PETSC ERROR: -nox [0]PETSC ERROR: -snes_monitor [0]PETSC ERROR: -ts_max_steps 2 [0]PETSC ERROR: -ts_type beuler [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint(a)mcs.anl.gov---------- application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 [unset]: aborting job: application called MPI_Abort(MPI_COMM_WORLD, 91) - process 0 I'm aware of Jed's commit e1244c6 which I thought was supposed to enable constant linear rhs without IFunction and IJacobian. Note that e1244c6 also added tests that apparently produced very different output (see e.g. src/ts/examples/tutorials/output/ex3_4.out). Is constant rhs without IFunction and IJacobian no longer supported or did it just breat at some point? Do I need to call TSRHSJacobianSetReuse in addition to the lines of code above? Is it better to treat such systems using TSComputeIFunctionLinear and TSComputeIJacobianConstant? Thanks, Dominic -- Dominic Meiser Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Karl Rupp 26 Feb '16

26 Feb '16

>> The alternative would be to use raw cuda pointers instead of cusp >> arrays for GPU memory in VecCUSP. That would be a fairly >> significant undertaking (certainly more than the 2-3 weeks Karli >> is estimating for getting the ViennaCL cuda backend in). > > Do you mean creating a new class VECCUDA in addition to VECCUSP and VECVIENNACL? This could be a solution for us. It would mean maybe refactoring MATAIJCUSPARSE to work with these new Vecs? I prefer to replace VECCUSP with e.g. VECCUDA (and eventually also rename VECCUSPARSE to VECCUDA to have a unified naming for all the things provided natively with the CUDA SDK) in order to reduce external dependencies. CUSP will provide matrices, preconditioners, etc. as before, but is only optional and thus less likely to cause installation troubles. Supporting VECCUSP and VECCUDA next to each other is going to be too much code duplication without any benefit. Even if we do provide VECCUDA, I still dislike the fact that we would have to maintain essentially the same code twice: One for CUDA, one for OpenCL/ViennaCL. With the ViennaCL bindings providing OpenMP and CUDA support soon, this also duplicates functionality. A possible 'fix' is to just use ViennaCL for CUDA+OpenCL+OpenMP and thus only maintain a single PETSc plugin for all three. However, I'm certainly too biased to be taken seriously here. > > If there is interest we can help in adding this stuff. What are your time constraints? Best regards, Karli

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Jose E. Roman 26 Feb '16

26 Feb '16

> El 25 feb 2016, a las 17:19, Dominic Meiser <dmeiser(a)txcorp.com> escribió: > > On Thu, Feb 25, 2016 at 01:13:01PM +0100, Jose E. Roman wrote: >> We are trying to do some GPU developments on the SLEPc side, and we would need a way of placing the array of a VECCUSP vector, providing the GPU address. Specifically, what we want to do is have a large Vec on GPU and slice it in several smaller Vecs. >> >> For the GetArray/RestoreArray we have all possibilities: >> - VecGetArray: gets the pointer to the buffer stored in CPU memory >> - VecCUSPGetArray*: returns a CUSPARRAY object that contains some info, including the buffer allocated in GPU memory >> - VecCUSPGetCUDAArray*: returns a raw pointer of the GPU buffer >> >> The problem comes with PlaceArray equivalents. Using VecPlaceArray we can provide a new pointer to CPU memory. We wanted to implement the equivalent thing for GPU, but we found difficulties due to Thrust. If we wanted to provide a VecCUSPPlaceCUDAArray the problem is that Thrust does not allow wrapping an exisiting GPU buffer with a CUSPARRAY object (when creating a CUSPARRAY it always allocates new memory). On the other hand, VecCUSPPlaceArray is possible to implement, but the problem is that one should provide a CUSPARRAY obtained from a VecCUSPGetArray* without modification (it is not possible to do pointer arithmetic with a CUSPARRAY). >> >> Any thoughts? >> > > I think your and Karli's analysis is correct, this is currently > not supported. Besides Karli's proposal to use ViennaCL's cuda > backend a different option might be to use cusp's array views. > These have a constructor for sub-ranges of other cusp arrays: > > https://github.com/cusplibrary/cusplibrary/blob/master/cusp/array1d.h#L409 > > However, enabling cusp array views in something like > VecCUSPPlaceArray is not immediately possible. The CUSPARRAY > type, which is currently hardwired to be > cusp::array1d<PetscScalar,cusp::device_memory>, would have to > become a template parameter. I'm not sure if we want to go down > that path. Yes, we do not like this. > > The alternative would be to use raw cuda pointers instead of cusp > arrays for GPU memory in VecCUSP. That would be a fairly > significant undertaking (certainly more than the 2-3 weeks Karli > is estimating for getting the ViennaCL cuda backend in). Do you mean creating a new class VECCUDA in addition to VECCUSP and VECVIENNACL? This could be a solution for us. It would mean maybe refactoring MATAIJCUSPARSE to work with these new Vecs? If there is interest we can help in adding this stuff. > > Cheers, > Dominic > > -- > Dominic Meiser > Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0

Re: [petsc-dev] Not possible to do a VecPlaceArray for veccusp
by Dominic Meiser 25 Feb '16

25 Feb '16

On Thu, Feb 25, 2016 at 01:13:01PM +0100, Jose E. Roman wrote: > We are trying to do some GPU developments on the SLEPc side, and we would need a way of placing the array of a VECCUSP vector, providing the GPU address. Specifically, what we want to do is have a large Vec on GPU and slice it in several smaller Vecs. > > For the GetArray/RestoreArray we have all possibilities: > - VecGetArray: gets the pointer to the buffer stored in CPU memory > - VecCUSPGetArray*: returns a CUSPARRAY object that contains some info, including the buffer allocated in GPU memory > - VecCUSPGetCUDAArray*: returns a raw pointer of the GPU buffer > > The problem comes with PlaceArray equivalents. Using VecPlaceArray we can provide a new pointer to CPU memory. We wanted to implement the equivalent thing for GPU, but we found difficulties due to Thrust. If we wanted to provide a VecCUSPPlaceCUDAArray the problem is that Thrust does not allow wrapping an exisiting GPU buffer with a CUSPARRAY object (when creating a CUSPARRAY it always allocates new memory). On the other hand, VecCUSPPlaceArray is possible to implement, but the problem is that one should provide a CUSPARRAY obtained from a VecCUSPGetArray* without modification (it is not possible to do pointer arithmetic with a CUSPARRAY). > > Any thoughts? > I think your and Karli's analysis is correct, this is currently not supported. Besides Karli's proposal to use ViennaCL's cuda backend a different option might be to use cusp's array views. These have a constructor for sub-ranges of other cusp arrays: https://github.com/cusplibrary/cusplibrary/blob/master/cusp/array1d.h#L409 However, enabling cusp array views in something like VecCUSPPlaceArray is not immediately possible. The CUSPARRAY type, which is currently hardwired to be cusp::array1d<PetscScalar,cusp::device_memory>, would have to become a template parameter. I'm not sure if we want to go down that path. The alternative would be to use raw cuda pointers instead of cusp arrays for GPU memory in VecCUSP. That would be a fairly significant undertaking (certainly more than the 2-3 weeks Karli is estimating for getting the ViennaCL cuda backend in). Cheers, Dominic -- Dominic Meiser Tech-X Corporation - 5621 Arapahoe Avenue - Boulder, CO 80303

1 0