[petsc-users] stdpar (nvfortran) managed memory + PETSc CUDA recommended interop pattern?
Junchao Zhang
junchao.zhang at gmail.com
Tue Jun 23 22:21:42 CDT 2026
Hi Edoardo,
Thanks for the example. First, I think we could easily add Fortran
bindings for VecCUDAGetArrayWrite() etc. We already have that for
VecGetArray().
We are sort of sloppy in logging every CPU-GPU transfer. We missed the
logging of a CpuToGpu transfer in MatSetValuesCOO. I fixed it in
branch jczhang/2026-06-23/fix-setvaluescoo-logging.
I also fixed the code and let petsc correctly detect the memory type of the
parameter coo_v[] in MatSetValuesCOO().
With this branch, I ran your test on a Linux machine with CUDA-12.8.
./ex01 -n 20000 -mat_type aijcusparse -log_view -log_view_gpu_time
coo_v pointer type as PETSc's COO path sees it:
coo_v 0x7f839a000000 : type=3 (Managed)
MatCreateVecs -> b type = seqcuda
||A*1||_2 = 1.41421
...
Event Count Time (sec) Flop
--- Global --- --- Stage ---- Total GPU - CpuToGpu - -
GpuToCpu - GPU
Max Ratio Max Ratio Max Ratio Mess AvgLen
Reduct %T %F %M %L %R %T %F %M %L %R Mflop/s Mflop/s Count Size
Count Size %F
---------------------------------------------------------------------------------------------------------------------------------------------------------------
--- Event Stage 0: Main Stage
MatMult 1 1.0 2.3852e-03 1.0 1.00e+05 1.0 0.0e+00 0.0e+00
0.0e+00 2 71 0 0 0 2 71 0 0 0 42 43 0 0.00e+00 0
0.00e+00 100
MatAssemblyBegin 1 1.0 1.4860e-06 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
MatAssemblyEnd 1 1.0 1.1236e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
*MatCUSPARSCopyTo 1 1.0 1.5919e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 1 0 0 0 0 1 0 0 0 0 0 0 1 8.00e-01 0
0.00e+00 0*
MatSetPreallCOO 1 1.0 3.0377e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 2 0 0 0 0 2 0 0 0 0 0 0 1 8.00e-01 0
0.00e+00 0
*MatSetValuesCOO 1 1.0 2.6547e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0*
cuBLAS Init 1 1.0 9.0025e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 7 0 0 0 0 7 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
DCtxCreate 2 1.0 4.4765e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
DCtxDestroy 2 1.0 1.6881e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
DCtxSetUp 2 1.0 1.1870e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
DCtxSetDevice 2 1.0 1.2694e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
VecNorm 1 1.0 1.0956e-03 1.0 4.00e+04 1.0 0.0e+00 0.0e+00
0.0e+00 1 29 0 0 0 1 29 0 0 0 37 38 0 0.00e+00 0
0.00e+00 100
VecSet 1 1.0 1.0189e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00
0.0e+00 8 0 0 0 0 8 0 0 0 0 0 0 0 0.00e+00 0
0.00e+00 0
You can see *MatSetValuesCOO* didn't incur CpuToGpu transfer. There was one
*MatCUSPARSCopyTo* since in MatSetPreallocationCOO(), we copied the
preallocated matrix data (zeros actually) on the host to the device.
Note that coo_v[] is not the matrix's data array on the device. petsc
MatAIJCUSPARSE matrix has its internal data storage on the device.
I think a really interesting questions is whether we can use
VecCreateSeqCUDAWithArrays(...,
managed_array, managed_array, &v) to create a CUDA vector sharing the same
CPU aray and GPU array, so that petsc users can seamlessly couple their
cuda managed data with petsc's VecCUDA.
We need to do more experiments.
--Junchao Zhang
On Mon, Jun 22, 2026 at 4:07 PM Edoardo alinovi <edoardo.alinovi at gmail.com>
wrote:
> Thanks Junchao, that clears it up.
>
> One open question: the CUDA device-array getters (VecCUDAGetArrayWrite,
> VecCreateMPICUDAWithArray) don't have Fortran stubs in my build. Is there a
> supported Fortran path (e.g. VecGetArrayAndMemType from Fortran), or are
> thin C shims the recommended approach?
>
> On the managed coo_v point: I rebuilt my PETSc with logging and ran the
> attached reproducer. With a cudaMallocManaged coo_v
> (cudaPointerGetAttributes → type = Managed) feeding MatSetValuesCOO on an
> AIJCUSPARSE matrix, -log_view shows MatSetValuesCOO with CpuToGpu Count =
> 0; the only H2D is the standard MatCUSPARSCopyTo (8 MB at n=2·10⁵). So I
> don't see an extra copy from the managed classification. The assembly is
> correct (MatMult checks out) and there's a single matrix to device copy.
>
> One caveat: I'm on WSL2, where cudaMemPrefetchAsync fails ("invalid device
> ordinal"), so coo_v stays host-resident. I couldn't exercise a
> device-resident managed coo_v, which is the case where a Device-only check
> might force a wasteful copy. Does that match your expectation, and is the
> device-resident-managed case worth me checking on a native-Linux box?
> Probably can double check next week on another machine.
>
> I am attaching a small test just in case it is useful to somebody out
> there.
>
> Build Info: PETSc 3.25.1, --with-cuda --with-cuda-arch=61, CUDA 12.6,
> gfortran-13.
>
> Thanks again,
>
> Edoardo
>
>
> Il giorno lun 22 giu 2026 alle ore 06:54 Junchao Zhang <
> junchao.zhang at gmail.com> ha scritto:
>
>> Hi Edoardo,
>> petsc/cuda backend doesn't use cuda managed memory, so I am just
>> trying to answer your questions based on my limited experience with it. See
>> the inlined answers.
>>
>>
>> On Sun, Jun 21, 2026 at 6:52 AM Edoardo alinovi <
>> edoardo.alinovi at gmail.com> wrote:
>>
>>> Hi PETSc friends!
>>>
>>> Hope you are all doing great.
>>>
>>> With my time I am porting my code to gpu. This is a really instructive
>>> project but it comes with several headaches as well.
>>>
>>> I have some valuable questions you might be able to answer that would
>>> help me a lot understanding how things work in the modern world:
>>>
>>> *My goal: *
>>> Assembly loops are *do concurrent* compiled with nvfortran -stdpar=gpu,
>>> so my field/work arrays live in CUDA managed memory (cudaMallocManaged).
>>> The linear solve uses PETSc with MATAIJCUSPARSE + (intended) VECCUDA. PETSc
>>> 3.25.1 (git), built with NVHPC 24.11 / CUDA 12.6; targets cc61 and cc86 (no
>>> HMM),. Iwant the whole outer iteration to stay device-resident.
>>>
>>> 1) Managed pointer as a VECCUDA device array. Can I wrap a
>>> cudaMallocManaged pointer as the device array of a Vec via
>>> VecCreateMPICUDAWithArray (or VecCUDAPlaceArray)? Will the offload-mask
>>> logic treat managed memory coherently, or assume a distinct host array and
>>> issue redundant H2D/D2H copies? Recommended pattern for stdpar-managed +
>>> PETSc-CUDA interop?
>>>
>>
>> Yes, I think you can use VecCreateMPICUDAWithArray() and friends.
>> petsc will just treat the array you provided as a device accessible array.
>> If you don't provide host arrays, petsc will malloc memory on the host to
>> mirror your cuda managed memory.
>>
>>
>>>
>>> 2) On-device RHS fill / solution read. Icurrently use
>>> VecSetValues/VecGetValues against host arrays, which flips the mask to CPU
>>> every iteration. With VECCUDA, is the recommended replacement
>>> VecCUDAGetArrayWrite (RHS) and VecCUDAGetArray (solution)? Offload-mask
>>> pitfalls?
>>>
>> I think you can do it. You just need to pretend your provided memory to
>> petsc is cude device memory. For example, you can call VecCUDAGetArrayWrite()
>> to get the array, then operate on host with the array, and then VecCUDAReturnArrayWrite().
>> You need to manage the synchronization yourself.
>>
>>
>>>
>>> 3) MatSetValuesCOO from device memory. Ipreallocate with
>>> MatSetPreallocationCOO and assemble with MatSetValuesCOO. For an
>>> AIJCUSPARSE matrix, may coo_v be a device/managed pointer? Memtype-aware
>>> contract, and from which version?
>>>
>>
>> petsc checks if the array coo_v is cudaMemoryTypeDevice. Perhaps we
>> should also check if it is cudaMemoryTypeManaged.
>>
>>
>>> 4) Vec type inference. If Iset only mat_type aijcusparse but leave the
>>> Vecs VECMPI, does KSP bounce the vectors each iteration, or promote them?
>>> Must I set vec_type cuda explicitly?
>>>
>> the vectors must also be device vectors. If you get the vectors from
>> MatGetVecs() or VecDuplicate() of the results, the vector types are
>> derived from matrix types (e.g., aijcusparse --> veccuda)
>>
>>
>>>
>>> Happy to share the configure line and a minimal reproducer.
>>>
>> Yes, please share one and issues you found, so we can better understand
>> the problem.
>>
>>
>>>
>>> Mega Thanks!
>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260623/7dbf44da/attachment.html>
More information about the petsc-users
mailing list