Let's see if we can lift this discussion up another level and "treat"
multi-core threading more specifically in the discussion (though Karl's subject
name is Unification approach for OpenMP/Threads/... he largely ignores the
multi-core/multi-socket aspect).
Abstractly a node has
1) a bunch of memories (some may be "nested" as caches "standing in" for parts
of larger caches which "stand in" for parts of "main memory". ) In general,
even without GPUs there are multiple memory sockets (though generally handled
by the OS as a single unified address space),
2) a bunch of compute "thingies". In general, even without GPUs there are
multiple CPUs, and each one of those likely has "regular" floating point units
plus SIMD units.
A) Shri has started coding up a runtime dispatch system for computations on
multiple cores (which hides differences between PThreads and OpenMP) that
(currently) assumes Vecs are stored in a single array (each thread accesses the
array pointer via VecGetArray() and then "its" part of the array by an offset.)
(BTW: what if each of these VecGetArray() triggered a copy up from a GPU,
probably a mess). When using PThreads Shir's model allows (to some degree) the
asynchronous launching of computational tasks.
B) We have a different dispatch system for using a single GPU accelerator via
CUDA that "automagically" handles copying data back and forth from memories via
VecXXXGetArray(). It is synchronous on the GetArray() in that is always blocks
on the GetArray() until the data is there and then moves on to the computation.
C) We are considering options for using OpenCL kernels.
D) We have not seriously considering utilizing both GPUs and core processors
for floating point intensive computations at the same time, either on the
"same" object computation or completely different object computations. (note
that DOE bought this huge machine at ORNL that seems to require this).
Ideally we'd have a "single" high performing programming model for utilizing
the resources of (1-2) regardless of details.
Now, lets go to Karl's "Part 1: Memory" which is a good place to start. In
PETSc we basically have two data types, a Vec which is relatively easy to
abstract about and a Mat which is not. Let's focus just on the Vec now because
Mat's are hard.
We need to "divide up" the computation on a Vec (or several Vecs and Mats)
so that the different compute "thingies" can work on their "piece", this
division of the computation naturally is associated with a "division" of the
data (the division may actually be only abstract with pthreads or it may be
concrete with two GPUs when "half" of the vector is copied to each GPU's memory
(sorry Jed, I agree with Karl that we likely shouldn't hide this issue behind
MPI)). The "division" is non-overlapping in simple cases (like axpy()) or may
require "ghosting" for sparse matrix-vector products (again the division my
only be abstract). With multi-memory-socket multi-core we actually divide the
vector data across physical memories but access it via virtual memory as not
divided up for ghost points etc. I think the "special cases" like virtual
memory make it harder for us to think about this abstractly then it should be.
In PETSc we use the abstract object IS to indicate parts of Vecs\footnote.
Thus if a computation requires part of a vector it is natural to pass into the
function the Vec AND THE IS indicating that part of the Vec needed. Note that
Shri's use of code such as i=trstarts[thread_id] is actually a particular type
of IS (hardwired for performance).
So, could we use a single kernel launcher for multi-core, CUDA, OpenCL based
on this principle? Then VecCUDAGetArray() type things would keep track of parts
of Vecs based on IS instead of all entries in the Vec. Similarly there would
be a VecMultiCoreGetArray(). Whenever possible the VecXXXGetArray() would not
require copies. As part of this model I'd also like to separate the "moving
needed data" part of the kernel from the "computation on the data" so that
everything doesn't block when data is being moved around.
Ok, how about moving this same model up to the MPI level? We already do this
with IS converted to VecScatter (for performance) for updating ghost points
(for matrix-vector products, for PDE ghost points etc) (note we can hide the
VecScatter inside the IS and have it created as needed).
Note I intend this to continue the conversation, not end it. Thoughts?
Barry
Footnote: Except when some people forget and make other unneeded complicated
constructs that reproduce the functionality of IS.
On Oct 6, 2012, at 9:09 AM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> On Sat, Oct 6, 2012 at 8:51 AM, Karl Rupp <rupp at mcs.anl.gov> wrote:
> Hmm, I thought that spptr is some 'special pointer' as commented in Mat, but
> not supposed to be a generic pointer to a derived class' datastructure (spptr
> is only injected with #define PETSC_HAVE_CUSP).
>
> Look at cuspvecimpl.h, for example.
>
> #undef __FUNCT__
> #define __FUNCT__ "VecCUSPGetArrayReadWrite"
> PETSC_STATIC_INLINE PetscErrorCode VecCUSPGetArrayReadWrite(Vec v,
> CUSPARRAY** a)
> {
> PetscErrorCode ierr;
>
> PetscFunctionBegin;
> *a = 0;
> ierr = VecCUSPCopyToGPU(v);CHKERRQ(ierr);
> *a = ((Vec_CUSP *)v->spptr)->GPUarray;
> PetscFunctionReturn(0);
> }
>
> Vec is following the convention from Mat where spptr points to the
> Mat_UMFPACK, Mat_SuperLU, Mat_SeqBSTRM, etc, which hold the extra information
> needed for that derived class (note that "derivation" for Mat is typically
> done at run time, after the object has been created, due to a call to
> MatGetFactor).
>
>
> With different devices, we could simply have a valid flag for each
> device. When someone does VecDevice1GetArrayWrite(), the flag for all
> other devices is marked invalid. When VecDevice2GetArrayRead() is
> called, the implementation copies from any valid device to device2.
> Packing all those flags as bits in a single int is perhaps convenient,
> but not necessary.
>
> I think that the most common way of handling GPUs will be an overlapping
> decomposition of the host array, similar to how a vector is distributed via
> MPI (locally owned, writeable, vs ghost values with read-only). Assigning the
> full vector exclusively to just one device is more a single-GPU scenario
> rather than a multi-GPU use case.
>
> Okay, the matrix will have to partition itself. What is the advantage of
> having a single CPU process addressing multiple GPUs? Why not use different
> MPI processes? (We can have the MPI processes sharing a node create a subcomm
> so they can decide which process is driving which device.)
>
>
> I think this stuff (which allows for segmenting the array on the device)
> can go in Vec_CUDA and Vec_OpenCL, basically just replacing the GPUarray
> member of Vec_CUSP. Why have a different PetscAcceleratorData struct?
>
> If spptr is intended to be a generic pointer to data of the derived class,
> then this is also a possiblity. However, this would lead to
> Vec_CUDA, Vec_OpenCL, and Vec_CUDA_OpenCL, with the number of implementations
> rapidly increasing as one may eventually add other frameworks. The
> PetscAcceleratorData would essentially allow for a unification of Vec_CUDA,
> Vec_OpenCL, and Vec_CUDA_OpenCL, avoiding code duplication problems.
>
> How would the user decide which device they wanted computation to run on?
> (Also, Is OpenCL really the right name in an environment where there may be
> multiple devices using OpenCL?) Currently, the type indicates where native
> operations should "prefer" to compute, copying data there when necessary. The
> Vec operations have different implementations for CUDA and OpenCL so I don't
> see the problem with making them different derived classes. If we wanted a
> hybrid CUDA/OpenCL class, it would contain the logic for deciding where to do
> things followed by dispatch into the device-specific implementation, thus it
> doesn't seem like duplication to me.
>
>
>
>
>
> Here, the PetscXYZHandleDescriptor holds
> - the memory handle,
> - the device ID the handles are valid for, and
> - a flag whether the data is valid
> (cf. valid_GPU_array, but with a much finer granularity).
> Additional metainformation such as index ranges can be extended as
> needed, cf. Vec_Seq vs Vec_MPI. Different types
> Petsc*HandleDescriptors are expected to be required because the
> various memory handle types are not guaranteed to have a particular
> maximum size among different accelerator platforms.
>
>
> It sounds like you want to support marking only part of an array as
> stale. We could could keep one top-level (_p_Vec) flag indicating
> whether the CPU part was current, then in the specific implementation
> (Vec_OpenCL), you can hold finer granularity. Then when
> vec->ops->UpdateCPUArray() is called, you can look at the finer
> granularity flags to copy only what needs to be copied.
>
>
> Yes, I also thought of such a top-level-flag. This is, however, rather an
> optimization flag (similar to what is done in VecGetArray for petscnative),
> so I refrained from a separate discussion.
>
> Aside from that, yes, I want to support parts of an array as stale, as the
> best multi-GPU use I've experienced so far is for block-based preconditioners
> (cf. Block-ILU-variants, parallel AMG flavors, etc.). A multi-GPU sparse
> matrix-vector product is handled similarly.
>
