Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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_sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 - 2022 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 #ifndef dealii_cuda_sparse_matrix_h
17 #define dealii_cuda_sparse_matrix_h
18 
19 #include <deal.II/base/config.h>
20 
22 
23 #include <iomanip>
24 
25 #ifdef DEAL_II_WITH_CUDA
26 # include <deal.II/base/cuda.h>
27 
28 # include <deal.II/lac/cuda_vector.h>
30 
31 # include <cusparse.h>
32 
34 
35 namespace CUDAWrappers
36 {
47  template <typename Number>
48  class SparseMatrix : public virtual Subscriptor
49  {
50  public:
54  using size_type = int;
55 
59  using value_type = Number;
60 
65  using real_type = Number;
66 
78  SparseMatrix();
79 
86  const ::SparseMatrix<Number> &sparse_matrix_host);
87 
93 
98 
102  ~SparseMatrix();
103 
107  SparseMatrix &
109 
113  SparseMatrix &
115 
121  void
123  const ::SparseMatrix<Number> &sparse_matrix_host);
134  size_type
135  m() const;
136 
141  size_type
142  n() const;
143 
149  std::size_t
150  n_nonzero_elements() const;
151 
162  template <class StreamType>
163  void
164  print(StreamType &out,
165  const bool across = false,
166  const bool diagonal_first = true) const;
167 
188  void
189  print_formatted(std::ostream & out,
190  const unsigned int precision = 3,
191  const bool scientific = true,
192  const unsigned int width = 0,
193  const char * zero_string = " ",
194  const double denominator = 1.) const;
204  SparseMatrix &
205  operator*=(const Number factor);
206 
210  SparseMatrix &
211  operator/=(const Number factor);
222  void
225 
231  void
234 
239  void
242 
248  void
251 
261  Number
264 
268  Number
272 
280  Number
296  Number
297  l1_norm() const;
298 
306  Number
307  linfty_norm() const;
308 
313  Number
314  frobenius_norm() const;
326  std::tuple<Number *, int *, int *, cusparseMatDescr_t, cusparseSpMatDescr_t>
327  get_cusparse_matrix() const;
330  private:
334  cusparseHandle_t cusparse_handle;
335 
339  int nnz;
340 
344  int n_rows;
345 
349  int n_cols;
350 
354  std::unique_ptr<Number[], void (*)(Number *)> val_dev;
355 
359  std::unique_ptr<int[], void (*)(int *)> column_index_dev;
360 
364  std::unique_ptr<int[], void (*)(int *)> row_ptr_dev;
365 
369  cusparseMatDescr_t descr;
370 
374  cusparseSpMatDescr_t sp_descr;
375  };
376 
377 
378 
379  template <typename Number>
380  inline typename SparseMatrix<Number>::size_type
382  {
383  return n_rows;
384  }
385 
386 
387 
388  template <typename Number>
389  inline typename SparseMatrix<Number>::size_type
391  {
392  return n_cols;
393  }
394 
395 
396 
397  template <typename Number>
398  inline std::size_t
400  {
401  return nnz;
402  }
403 
404 
405 
406  template <typename Number>
407  template <class StreamType>
408  inline void
410  const bool across,
411  const bool diagonal_first) const
412  {
413  Assert(column_index_dev != nullptr, ExcNotInitialized());
414  Assert(val_dev != nullptr, ExcNotInitialized());
415  Assert(row_ptr_dev != nullptr, ExcNotInitialized());
416 
417  std::vector<int> rows(n_rows + 1);
418  std::vector<int> cols(nnz);
419  std::vector<double> val(nnz);
420  Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
421  Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
422  Utilities::CUDA::copy_to_host(val_dev.get(), val);
423 
424  bool has_diagonal = false;
425  Number diagonal = Number();
426 
427  for (size_type i = 0; i < n_rows; ++i)
428  {
429  if (diagonal_first)
430  {
431  // find the diagonal and print if it exists
432  for (size_type j = rows[i]; j < rows[i + 1] && cols[j] <= i; ++j)
433  {
434  if (i == cols[j])
435  {
436  diagonal = val[j];
437  has_diagonal = true;
438  if (across)
439  out << ' ' << i << ',' << i << ':' << diagonal;
440  else
441  out << '(' << i << ',' << i << ") " << diagonal
442  << std::endl;
443  break;
444  }
445  }
446  }
447  for (size_type j = rows[i]; j < rows[i + 1]; ++j)
448  {
449  if (has_diagonal && i == cols[j])
450  continue;
451  if (across)
452  out << ' ' << i << ',' << cols[j] << ':' << val[j];
453  else
454  out << '(' << i << ',' << cols[j] << ") " << val[j] << std::endl;
455  }
456  }
457  if (across)
458  out << std::endl;
459  }
460 
461 
462 
463  template <typename Number>
464  void
466  const unsigned int precision,
467  const bool scientific,
468  const unsigned int width_,
469  const char * zero_string,
470  const double denominator) const
471  {
472  Assert(column_index_dev != nullptr, ExcNotInitialized());
473  Assert(val_dev != nullptr, ExcNotInitialized());
474  Assert(row_ptr_dev != nullptr, ExcNotInitialized());
475 
476  std::vector<int> rows(n_rows + 1);
477  std::vector<int> cols(nnz);
478  std::vector<Number> val(nnz);
479  Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
480  Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
481  Utilities::CUDA::copy_to_host(val_dev.get(), val);
482 
483  unsigned int width = width_;
484 
485  std::ios::fmtflags old_flags = out.flags();
486  unsigned int old_precision = out.precision(precision);
487 
488  if (scientific)
489  {
490  out.setf(std::ios::scientific, std::ios::floatfield);
491  if (!width)
492  width = precision + 7;
493  }
494  else
495  {
496  out.setf(std::ios::fixed, std::ios::floatfield);
497  if (!width)
498  width = precision + 2;
499  }
500 
501  for (size_type i = 0; i < n_rows; ++i)
502  {
503  size_type j = rows[i];
504  for (size_type k = 0; k < n_cols; ++k)
505  {
506  if (k == cols[j])
507  {
508  out << std::setw(width) << val[j] * Number(denominator) << ' ';
509  ++j;
510  }
511  else
512  out << std::setw(width) << zero_string << ' ';
513  }
514  out << std::endl;
515  };
516  AssertThrow(out.fail() == false, ExcIO());
517 
518  // reset output format
519  out.precision(old_precision);
520  out.flags(old_flags);
521  }
522 } // namespace CUDAWrappers
523 
525 
526 #endif
527 #endif
void vmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
std::tuple< Number *, int *, int *, cusparseMatDescr_t, cusparseSpMatDescr_t > get_cusparse_matrix() const
void Tvmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
std::unique_ptr< int[], void(*)(int *)> column_index_dev
void Tvmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
std::unique_ptr< int[], void(*)(int *)> row_ptr_dev
void reinit(Utilities::CUDA::Handle &handle, const ::SparseMatrix< Number > &sparse_matrix_host)
Number matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
SparseMatrix & operator*=(const Number factor)
cusparseSpMatDescr_t sp_descr
Number residual(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &x, const LinearAlgebra::CUDAWrappers::Vector< Number > &b) const
SparseMatrix & operator/=(const Number factor)
SparseMatrix & operator=(CUDAWrappers::SparseMatrix< Number > &&)
std::unique_ptr< Number[], void(*)(Number *)> val_dev
void vmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
SparseMatrix(const CUDAWrappers::SparseMatrix< Number > &)=delete
SparseMatrix & operator=(const CUDAWrappers::SparseMatrix< Number > &)=delete
std::size_t n_nonzero_elements() const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
Number matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector< Number > &u, const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcIO()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
@ diagonal
Matrix is diagonal.
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition: cuda.h:132