Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
cuda_sparse_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2018 - 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#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
30
31# include <cusparse.h>
32
34
35namespace 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
79
86 const ::SparseMatrix<Number> &sparse_matrix_host);
87
93
98
103
109
115
121 void
123 const ::SparseMatrix<Number> &sparse_matrix_host);
135 m() const;
136
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;
205 operator*=(const Number factor);
206
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
345
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
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
SparseMatrix & operator=(const CUDAWrappers::SparseMatrix< Number > &)=delete
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition cuda.h:132