Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
cuda_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2018 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 
21 #include <deal.II/base/subscriptor.h>
22 
23 #include <iomanip>
24 
25 #ifdef DEAL_II_COMPILER_CUDA_AWARE
26 # include <deal.II/base/cuda.h>
27 
28 # include <deal.II/lac/cuda_vector.h>
29 # include <deal.II/lac/sparse_matrix.h>
30 
31 # include <cusparse.h>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 namespace CUDAWrappers
36 {
49  template <typename Number>
50  class SparseMatrix : public virtual Subscriptor
51  {
52  public:
56  using size_type = int;
57 
61  using value_type = Number;
62 
67  using real_type = Number;
68 
80  SparseMatrix();
81 
88  const ::SparseMatrix<Number> &sparse_matrix_host);
89 
95 
100 
104  ~SparseMatrix();
105 
109  SparseMatrix &
111 
115  SparseMatrix &
117 
123  void
125  const ::SparseMatrix<Number> &sparse_matrix_host);
127 
136  size_type
137  m() const;
138 
143  size_type
144  n() const;
145 
151  std::size_t
152  n_nonzero_elements() const;
153 
164  template <class StreamType>
165  void
166  print(StreamType &out,
167  const bool across = false,
168  const bool diagonal_first = true) const;
169 
190  void
191  print_formatted(std::ostream & out,
192  const unsigned int precision = 3,
193  const bool scientific = true,
194  const unsigned int width = 0,
195  const char * zero_string = " ",
196  const double denominator = 1.) const;
198 
206  SparseMatrix &
207  operator*=(const Number factor);
208 
212  SparseMatrix &
213  operator/=(const Number factor);
215 
224  void
227 
233  void
236 
241  void
244 
250  void
253 
263  Number
266 
270  Number
274 
282  Number
287 
298  Number
299  l1_norm() const;
300 
308  Number
309  linfty_norm() const;
310 
315  Number
316  frobenius_norm() const;
318 
328  std::tuple<Number *, int *, int *, cusparseMatDescr_t>
329  get_cusparse_matrix() const;
331 
332  private:
336  cusparseHandle_t cusparse_handle;
337 
341  int nnz;
342 
346  int n_rows;
347 
351  int n_cols;
352 
356  std::unique_ptr<Number[], void (*)(Number *)> val_dev;
357 
361  std::unique_ptr<int[], void (*)(int *)> column_index_dev;
362 
366  std::unique_ptr<int[], void (*)(int *)> row_ptr_dev;
367 
371  cusparseMatDescr_t descr;
372  };
373 
374 
375 
376  template <typename Number>
379  {
380  return n_rows;
381  }
382 
383 
384 
385  template <typename Number>
388  {
389  return n_cols;
390  }
391 
392 
393 
394  template <typename Number>
395  inline std::size_t
397  {
398  return nnz;
399  }
400 
401 
402 
403  template <typename Number>
404  template <class StreamType>
405  inline void
407  const bool across,
408  const bool diagonal_first) const
409  {
410  Assert(column_index_dev != nullptr, ExcNotInitialized());
411  Assert(val_dev != nullptr, ExcNotInitialized());
412  Assert(row_ptr_dev != nullptr, ExcNotInitialized());
413 
414  std::vector<int> rows(n_rows + 1);
415  std::vector<int> cols(nnz);
416  std::vector<double> val(nnz);
417  Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
418  Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
419  Utilities::CUDA::copy_to_host(val_dev.get(), val);
420 
421  bool has_diagonal = false;
422  Number diagonal = Number();
423 
424  for (size_type i = 0; i < n_rows; ++i)
425  {
426  if (diagonal_first)
427  {
428  // find the diagonal and print if it exists
429  for (size_type j = rows[i]; j < rows[i + 1] && cols[j] <= i; ++j)
430  {
431  if (i == cols[j])
432  {
433  diagonal = val[j];
434  has_diagonal = true;
435  if (across)
436  out << ' ' << i << ',' << i << ':' << diagonal;
437  else
438  out << '(' << i << ',' << i << ") " << diagonal
439  << std::endl;
440  break;
441  }
442  }
443  }
444  for (size_type j = rows[i]; j < rows[i + 1]; ++j)
445  {
446  if (has_diagonal && i == cols[j])
447  continue;
448  if (across)
449  out << ' ' << i << ',' << cols[j] << ':' << val[j];
450  else
451  out << "(" << i << "," << cols[j] << ") " << val[j] << std::endl;
452  }
453  }
454  if (across)
455  out << std::endl;
456  }
457 
458 
459 
460  template <typename Number>
461  void
463  const unsigned int precision,
464  const bool scientific,
465  const unsigned int width_,
466  const char * zero_string,
467  const double denominator) const
468  {
469  Assert(column_index_dev != nullptr, ExcNotInitialized());
470  Assert(val_dev != nullptr, ExcNotInitialized());
471  Assert(row_ptr_dev != nullptr, ExcNotInitialized());
472 
473  std::vector<int> rows(n_rows + 1);
474  std::vector<int> cols(nnz);
475  std::vector<Number> val(nnz);
476  Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
477  Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
478  Utilities::CUDA::copy_to_host(val_dev.get(), val);
479 
480  unsigned int width = width_;
481 
482  std::ios::fmtflags old_flags = out.flags();
483  unsigned int old_precision = out.precision(precision);
484 
485  if (scientific)
486  {
487  out.setf(std::ios::scientific, std::ios::floatfield);
488  if (!width)
489  width = precision + 7;
490  }
491  else
492  {
493  out.setf(std::ios::fixed, std::ios::floatfield);
494  if (!width)
495  width = precision + 2;
496  }
497 
498  for (size_type i = 0; i < n_rows; ++i)
499  {
500  size_type j = rows[i];
501  for (size_type k = 0; k < n_cols; ++k)
502  {
503  if (k == cols[j])
504  {
505  out << std::setw(width) << val[j] * Number(denominator) << ' ';
506  ++j;
507  }
508  else
509  out << std::setw(width) << zero_string << ' ';
510  }
511  out << std::endl;
512  };
513  AssertThrow(out, ExcIO());
514 
515  // reset output format
516  out.precision(old_precision);
517  out.flags(old_flags);
518  }
519 } // namespace CUDAWrappers
520 
521 DEAL_II_NAMESPACE_CLOSE
522 
523 #endif
524 #endif
Number matrix_norm_square(const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
static ::ExceptionBase & ExcIO()
std::tuple< Number *, int *, int *, cusparseMatDescr_t > get_cusparse_matrix() const
Number matrix_scalar_product(const LinearAlgebra::CUDAWrappers::Vector< Number > &u, const LinearAlgebra::CUDAWrappers::Vector< Number > &v) const
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
SparseMatrix & operator*=(const Number factor)
Number frobenius_norm() const
SparseMatrix & operator/=(const Number factor)
std::unique_ptr< Number[], void(*)(Number *)> val_dev
std::unique_ptr< int[], void(*)(int *)> column_index_dev
Number residual(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &x, const LinearAlgebra::CUDAWrappers::Vector< Number > &b) const
void reinit(Utilities::CUDA::Handle &handle, const ::SparseMatrix< Number > &sparse_matrix_host)
Number linfty_norm() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::unique_ptr< int[], void(*)(int *)> row_ptr_dev
SparseMatrix & operator=(CUDAWrappers::SparseMatrix< Number > &&)
void Tvmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
void copy_to_host(const T *pointer_dev, std::vector< T > &vector_host)
Definition: cuda.h:131
void vmult_add(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) const
void Tvmult(LinearAlgebra::CUDAWrappers::Vector< Number > &dst, const LinearAlgebra::CUDAWrappers::Vector< Number > &src) 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
void vmult(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::size_t n_nonzero_elements() const