Reference documentation for deal.II version GIT f6a5d312c9 2023-10-04 08:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
cuda_vector.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/cuda.h>
17 #include <deal.II/base/cuda_size.h>
19 
24 
25 #include <cmath>
26 
27 #ifdef DEAL_II_WITH_CUDA
28 
30 
31 namespace LinearAlgebra
32 {
33  namespace CUDAWrappers
34  {
37 
38 
39 
40  template <typename Number>
42  : val(nullptr, Utilities::CUDA::delete_device_data<Number>)
43  , n_elements(0)
44  {}
45 
46 
47 
48  template <typename Number>
50  : val(Utilities::CUDA::allocate_device_data<Number>(V.n_elements),
52  , n_elements(V.n_elements)
53  {
54  // Copy the values.
55  const cudaError_t error_code = cudaMemcpy(val.get(),
56  V.val.get(),
57  n_elements * sizeof(Number),
58  cudaMemcpyDeviceToDevice);
59  AssertCuda(error_code);
60  }
61 
62 
63 
64  template <typename Number>
67  {
68  if (n_elements < V.n_elements)
69  reinit(V.n_elements, true);
70  else
71  n_elements = V.n_elements;
72 
73  // Copy the values.
74  const cudaError_t error_code = cudaMemcpy(val.get(),
75  V.val.get(),
76  n_elements * sizeof(Number),
77  cudaMemcpyDeviceToDevice);
78  AssertCuda(error_code);
79 
80  return *this;
81  }
82 
83 
84 
85  template <typename Number>
87  : val(nullptr, Utilities::CUDA::delete_device_data<Number>)
88  , n_elements(0)
89  {
90  reinit(n, false);
91  }
92 
93 
94 
95  template <typename Number>
96  void
97  Vector<Number>::reinit(const size_type n, const bool omit_zeroing_entries)
98  {
99  // Resize the underlying array if necessary
100  if (n == 0)
101  val.reset();
102  else if (n != n_elements)
103  val.reset(Utilities::CUDA::allocate_device_data<Number>(n));
104 
105  // If necessary set the elements to zero
106  if (omit_zeroing_entries == false)
107  {
108  const cudaError_t error_code =
109  cudaMemset(val.get(), 0, n * sizeof(Number));
110  AssertCuda(error_code);
111  }
112  n_elements = n;
113  }
114 
115 
116 
117  template <typename Number>
118  void
120  const bool omit_zeroing_entries)
121  {
122  reinit(V.size(), omit_zeroing_entries);
123  }
124 
125 
126 
127  template <typename Number>
128  void
130  const ReadWriteVector<Number> &V,
131  const VectorOperation::values operation,
132  const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase> &)
133  {
134  if (operation == VectorOperation::insert)
135  {
136  const cudaError_t error_code = cudaMemcpy(val.get(),
137  V.begin(),
138  n_elements * sizeof(Number),
139  cudaMemcpyHostToDevice);
140  AssertCuda(error_code);
141  }
142  else if (operation == VectorOperation::add)
143  {
144  // Create a temporary vector on the device
145  Number *tmp;
146  cudaError_t error_code =
147  cudaMalloc(&tmp, n_elements * sizeof(Number));
148  AssertCuda(error_code);
149 
150  // Copy the vector from the host to the temporary vector on the device
151  error_code = cudaMemcpy(tmp,
152  V.begin(),
153  n_elements * sizeof(Number),
154  cudaMemcpyHostToDevice);
155  AssertCuda(error_code);
156 
157  // Add the two vectors
158  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
159 
160  kernel::vector_bin_op<Number, kernel::Binop_Addition>
161  <<<n_blocks, block_size>>>(val.get(), tmp, n_elements);
163 
164  // Delete the temporary vector
166  }
167  else
168  AssertThrow(false, ExcNotImplemented());
169  }
170 
171 
172 
173  template <typename Number>
176  {
177  Assert(s == Number(), ExcMessage("Only 0 can be assigned to a vector."));
178  (void)s;
179 
180  const cudaError_t error_code =
181  cudaMemset(val.get(), 0, n_elements * sizeof(Number));
182  AssertCuda(error_code);
183 
184  return *this;
185  }
186 
187 
188 
189  template <typename Number>
191  Vector<Number>::operator*=(const Number factor)
192  {
193  AssertIsFinite(factor);
194  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
195  kernel::vec_scale<Number>
196  <<<n_blocks, block_size>>>(val.get(), factor, n_elements);
198 
199  return *this;
200  }
201 
202 
203 
204  template <typename Number>
206  Vector<Number>::operator/=(const Number factor)
207  {
208  AssertIsFinite(factor);
209  Assert(factor != Number(0.), ExcZero());
210  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
211  kernel::vec_scale<Number>
212  <<<n_blocks, block_size>>>(val.get(), 1. / factor, n_elements);
214 
215  return *this;
216  }
217 
218 
219 
220  template <typename Number>
223  {
224  Assert(V.size() == this->size(),
225  ExcMessage(
226  "Cannot add two vectors with different numbers of elements"));
227 
228  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
229 
230  kernel::vector_bin_op<Number, kernel::Binop_Addition>
231  <<<n_blocks, block_size>>>(val.get(), V.val.get(), n_elements);
233 
234  return *this;
235  }
236 
237 
238 
239  template <typename Number>
242  {
243  Assert(V.size() == this->size(),
244  ExcMessage(
245  "Cannot add two vectors with different numbers of elements."));
246 
247  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
248 
249  kernel::vector_bin_op<Number, kernel::Binop_Subtraction>
250  <<<n_blocks, block_size>>>(val.get(), V.val.get(), n_elements);
252 
253  return *this;
254  }
255 
256 
257 
258  template <typename Number>
259  Number
261  {
262  Assert(V.size() == this->size(),
263  ExcMessage(
264  "Cannot add two vectors with different numbers of elements"));
265 
266  Number *result_device;
267  cudaError_t error_code =
268  cudaMalloc(&result_device, n_elements * sizeof(Number));
269  AssertCuda(error_code);
270  error_code = cudaMemset(result_device, 0, sizeof(Number));
271 
272  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
273  kernel::double_vector_reduction<Number, kernel::DotProduct<Number>>
274  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
275  val.get(),
276  V.val.get(),
277  static_cast<unsigned int>(
278  n_elements));
279 
280  // Copy the result back to the host
281  Number result;
282  error_code = cudaMemcpy(&result,
283  result_device,
284  sizeof(Number),
285  cudaMemcpyDeviceToHost);
286  AssertCuda(error_code);
287  // Free the memory on the device
288  Utilities::CUDA::free(result_device);
289 
290  return result;
291  }
292 
293 
294 
295  template <typename Number>
296  void
297  Vector<Number>::add(const Number a)
298  {
299  AssertIsFinite(a);
300  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
301  kernel::vec_add<Number>
302  <<<n_blocks, block_size>>>(val.get(), a, n_elements);
304  }
305 
306 
307 
308  template <typename Number>
309  void
310  Vector<Number>::add(const Number a, const Vector<Number> &V)
311  {
312  AssertIsFinite(a);
313 
314  Assert(V.size() == this->size(),
315  ExcMessage(
316  "Cannot add two vectors with different numbers of elements."));
317 
318  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
319  kernel::add_aV<Number><<<dim3(n_blocks, 1), dim3(block_size)>>>(
320  val.get(), a, V.val.get(), n_elements);
322  }
323 
324 
325 
326  template <typename Number>
327  void
328  Vector<Number>::add(const Number a,
329  const Vector<Number> &V,
330  const Number b,
331  const Vector<Number> &W)
332  {
333  AssertIsFinite(a);
334  AssertIsFinite(b);
335 
336  Assert(V.size() == this->size(),
337  ExcMessage(
338  "Cannot add two vectors with different numbers of elements."));
339 
340  Assert(W.size() == this->size(),
341  ExcMessage(
342  "Cannot add two vectors with different numbers of elements."));
343 
344  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
345  kernel::add_aVbW<Number><<<dim3(n_blocks, 1), dim3(block_size)>>>(
346  val.get(), a, V.val.get(), b, W.val.get(), n_elements);
348  }
349 
350 
351 
352  template <typename Number>
353  void
354  Vector<Number>::sadd(const Number s,
355  const Number a,
356  const Vector<Number> &V)
357  {
358  AssertIsFinite(s);
359  AssertIsFinite(a);
360 
361  Assert(V.size() == this->size(),
362  ExcMessage(
363  "Cannot add two vectors with different numbers of elements."));
364 
365  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
366  kernel::sadd<Number><<<dim3(n_blocks, 1), dim3(block_size)>>>(
367  s, val.get(), a, V.val.get(), n_elements);
369  }
370 
371 
372 
373  template <typename Number>
374  void
375  Vector<Number>::scale(const Vector<Number> &scaling_factors)
376  {
377  Assert(scaling_factors.size() == this->size(),
378  ExcMessage(
379  "Cannot scale two vectors with different numbers of elements."));
380 
381  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
382  kernel::scale<Number>
383  <<<dim3(n_blocks, 1), dim3(block_size)>>>(val.get(),
384  scaling_factors.val.get(),
385  n_elements);
387  }
388 
389 
390 
391  template <typename Number>
392  void
393  Vector<Number>::equ(const Number a, const Vector<Number> &V)
394  {
395  AssertIsFinite(a);
396 
397  Assert(
398  V.size() == this->size(),
399  ExcMessage(
400  "Cannot assign two vectors with different numbers of elements."));
401 
402  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
403  kernel::equ<Number><<<dim3(n_blocks, 1), dim3(block_size)>>>(val.get(),
404  a,
405  V.val.get(),
406  n_elements);
408  }
409 
410 
411 
412  template <typename Number>
413  bool
415  {
416  return (linfty_norm() == 0) ? true : false;
417  }
418 
419 
420 
421  template <typename Number>
424  {
425  Number *result_device;
426  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
427  AssertCuda(error_code);
428  error_code = cudaMemset(result_device, 0, sizeof(Number));
429 
430  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
431  kernel::reduction<Number, kernel::ElemSum<Number>>
432  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
433  val.get(),
434  n_elements);
435 
436  // Copy the result back to the host
437  Number result;
438  error_code = cudaMemcpy(&result,
439  result_device,
440  sizeof(Number),
441  cudaMemcpyDeviceToHost);
442  AssertCuda(error_code);
443  // Free the memory on the device
444  Utilities::CUDA::free(result_device);
445 
446  return result /
447  static_cast<typename Vector<Number>::value_type>(n_elements);
448  }
449 
450 
451 
452  template <typename Number>
455  {
456  Number *result_device;
457  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
458  AssertCuda(error_code);
459  error_code = cudaMemset(result_device, 0, sizeof(Number));
460 
461  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
462  kernel::reduction<Number, kernel::L1Norm<Number>>
463  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
464  val.get(),
465  n_elements);
466 
467  // Copy the result back to the host
468  Number result;
469  error_code = cudaMemcpy(&result,
470  result_device,
471  sizeof(Number),
472  cudaMemcpyDeviceToHost);
473  AssertCuda(error_code);
474  // Free the memory on the device
475  Utilities::CUDA::free(result_device);
476 
477  return result;
478  }
479 
480 
481 
482  template <typename Number>
485  {
486  return std::sqrt(norm_sqr());
487  }
488 
489 
490 
491  template <typename Number>
494  {
495  return (*this) * (*this);
496  }
497 
498 
499 
500  template <typename Number>
503  {
504  Number *result_device;
505  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
506  AssertCuda(error_code);
507  error_code = cudaMemset(result_device, 0, sizeof(Number));
508 
509  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
510  kernel::reduction<Number, kernel::LInfty<Number>>
511  <<<dim3(n_blocks, 1), dim3(block_size)>>>(result_device,
512  val.get(),
513  n_elements);
514 
515  // Copy the result back to the host
516  Number result;
517  error_code = cudaMemcpy(&result,
518  result_device,
519  sizeof(Number),
520  cudaMemcpyDeviceToHost);
521  AssertCuda(error_code);
522  // Free the memory on the device
523  Utilities::CUDA::free(result_device);
524 
525  return result;
526  }
527 
528 
529 
530  template <typename Number>
531  Number
533  const Vector<Number> &V,
534  const Vector<Number> &W)
535  {
536  AssertIsFinite(a);
537 
538  Assert(V.size() == this->size(),
539  ExcMessage("Vector V has the wrong size."));
540  Assert(W.size() == this->size(),
541  ExcMessage("Vector W has the wrong size."));
542 
543  Number *result_device;
544  cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
545  AssertCuda(error_code);
546  error_code = cudaMemset(result_device, 0, sizeof(Number));
547  AssertCuda(error_code);
548 
549  const int n_blocks = 1 + (n_elements - 1) / (chunk_size * block_size);
550  kernel::add_and_dot<Number><<<dim3(n_blocks, 1), dim3(block_size)>>>(
551  result_device, val.get(), V.val.get(), W.val.get(), a, n_elements);
552 
553  Number result;
554  error_code = cudaMemcpy(&result,
555  result_device,
556  sizeof(Number),
557  cudaMemcpyDeviceToHost);
558  Utilities::CUDA::free(result_device);
559 
560  return result;
561  }
562 
563 
564 
565  template <typename Number>
566  void
567  Vector<Number>::print(std::ostream &out,
568  const unsigned int precision,
569  const bool scientific,
570  const bool) const
571  {
572  AssertThrow(out.fail() == false, ExcIO());
573  std::ios::fmtflags old_flags = out.flags();
574  unsigned int old_precision = out.precision(precision);
575 
576  out.precision(precision);
577  if (scientific)
578  out.setf(std::ios::scientific, std::ios::floatfield);
579  else
580  out.setf(std::ios::fixed, std::ios::floatfield);
581 
582  out << "IndexSet: ";
583  complete_index_set(n_elements).print(out);
584  out << std::endl;
585 
586  // Copy the vector to the host
587  std::vector<Number> cpu_val(n_elements);
588  Utilities::CUDA::copy_to_host(val.get(), cpu_val);
589  for (unsigned int i = 0; i < n_elements; ++i)
590  out << cpu_val[i] << std::endl;
591  out << std::flush;
592 
593  AssertThrow(out.fail() == false, ExcIO());
594  // reset output format
595  out.flags(old_flags);
596  out.precision(old_precision);
597  }
598 
599 
600 
601  template <typename Number>
602  std::size_t
604  {
605  std::size_t memory = sizeof(*this);
606  memory += sizeof(Number) * static_cast<std::size_t>(n_elements);
607 
608  return memory;
609  }
610 
611 
612 
613  // Explicit Instanationation
614  template class Vector<float>;
615  template class Vector<double>;
616  } // namespace CUDAWrappers
617 } // namespace LinearAlgebra
618 
620 
621 #endif
void print(StreamType &out) const
Definition: index_set.h:2009
void equ(const Number a, const Vector< Number > &V)
Definition: cuda_vector.cc:393
types::global_dof_index size_type
Definition: cuda_vector.h:60
void sadd(const Number s, const Number a, const Vector< Number > &V)
Definition: cuda_vector.cc:354
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
Definition: cuda_vector.cc:532
std::unique_ptr< Number[], void(*)(Number *)> val
Definition: cuda_vector.h:347
typename numbers::NumberTraits< Number >::real_type real_type
Definition: cuda_vector.h:61
Vector< Number > & operator*=(const Number factor)
Definition: cuda_vector.cc:191
void import_elements(const ReadWriteVector< Number > &V, const VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern={})
Definition: cuda_vector.cc:129
Vector & operator=(const Vector< Number > &v)
Definition: cuda_vector.cc:66
Vector< Number > & operator+=(const Vector< Number > &V)
Definition: cuda_vector.cc:222
void scale(const Vector< Number > &scaling_factors)
Definition: cuda_vector.cc:375
void print(std::ostream &out, const unsigned int precision=2, const bool scientific=true, const bool across=true) const
Definition: cuda_vector.cc:567
Number operator*(const Vector< Number > &V) const
Definition: cuda_vector.cc:260
void reinit(const size_type n, const bool omit_zeroing_entries=false)
Definition: cuda_vector.cc:97
Vector< Number > & operator-=(const Vector< Number > &V)
Definition: cuda_vector.cc:241
std::size_t memory_consumption() const
Definition: cuda_vector.cc:603
Vector< Number > & operator/=(const Number factor)
Definition: cuda_vector.cc:206
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcZero()
#define AssertCudaKernel()
Definition: exceptions.h:1997
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define AssertCuda(error_code)
Definition: exceptions.h:1942
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1124
constexpr int chunk_size
Definition: cuda_size.h:35
constexpr int block_size
Definition: cuda_size.h:29
static const char V
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
Default CUDA
Definition: memory_space.h:53
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number * allocate_device_data(const std::size_t size)
Definition: cuda.h:109
void delete_device_data(Number *device_ptr) noexcept
Definition: cuda.h:121
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition: cuda.h:132
void free(T *&pointer)
Definition: cuda.h:97
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3061