Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00: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\}}\)
Loading...
Searching...
No Matches
cuda_sparse_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2018 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_cuda_sparse_matrix_h
16#define dealii_cuda_sparse_matrix_h
17
18#include <deal.II/base/config.h>
19
21
22#include <iomanip>
23
24#ifdef DEAL_II_WITH_CUDA
25# include <deal.II/base/cuda.h>
26
29
30# include <cusparse.h>
31
33
34namespace CUDAWrappers
35{
46 template <typename Number>
47 class SparseMatrix : public virtual Subscriptor
48 {
49 public:
53 using size_type = int;
54
58 using value_type = Number;
59
64 using real_type = Number;
65
78
85 const ::SparseMatrix<Number> &sparse_matrix_host);
86
92
97
102
108
114
120 void
122 const ::SparseMatrix<Number> &sparse_matrix_host);
134 m() const;
135
141 n() const;
142
148 std::size_t
149 n_nonzero_elements() const;
150
161 template <typename StreamType>
162 void
163 print(StreamType &out,
164 const bool across = false,
165 const bool diagonal_first = true) const;
166
189 void
190 print_formatted(std::ostream &out,
191 const unsigned int precision = 3,
192 const bool scientific = true,
193 const unsigned int width = 0,
194 const char *zero_string = " ",
195 const double denominator = 1.,
196 const char *separator = " ") const;
207 operator*=(const Number factor);
208
213 operator/=(const Number factor);
224 void
227
233 void
236
241 void
244
250 void
253
263 Number
266
270 Number
274
282 Number
298 Number
299 l1_norm() const;
300
308 Number
309 linfty_norm() const;
310
315 Number
316 frobenius_norm() const;
328 std::tuple<Number *, int *, int *, cusparseMatDescr_t, cusparseSpMatDescr_t>
329 get_cusparse_matrix() const;
332 private:
336 cusparseHandle_t cusparse_handle;
337
341 int nnz;
342
347
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
376 cusparseSpMatDescr_t sp_descr;
377 };
378
379
380
381 template <typename Number>
382 inline typename SparseMatrix<Number>::size_type
384 {
385 return n_rows;
386 }
387
388
389
390 template <typename Number>
391 inline typename SparseMatrix<Number>::size_type
393 {
394 return n_cols;
395 }
396
397
398
399 template <typename Number>
400 inline std::size_t
402 {
403 return nnz;
404 }
405
406
407
408 template <typename Number>
409 template <typename StreamType>
410 inline void
412 const bool across,
413 const bool diagonal_first) const
414 {
415 Assert(column_index_dev != nullptr, ExcNotInitialized());
416 Assert(val_dev != nullptr, ExcNotInitialized());
417 Assert(row_ptr_dev != nullptr, ExcNotInitialized());
418
419 std::vector<int> rows(n_rows + 1);
420 std::vector<int> cols(nnz);
421 std::vector<double> val(nnz);
422 Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
423 Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
424 Utilities::CUDA::copy_to_host(val_dev.get(), val);
425
426 bool has_diagonal = false;
427 Number diagonal = Number();
428
429 for (size_type i = 0; i < n_rows; ++i)
430 {
431 if (diagonal_first)
432 {
433 // find the diagonal and print if it exists
434 for (size_type j = rows[i]; j < rows[i + 1] && cols[j] <= i; ++j)
435 {
436 if (i == cols[j])
437 {
438 diagonal = val[j];
439 has_diagonal = true;
440 if (across)
441 out << ' ' << i << ',' << i << ':' << diagonal;
442 else
443 out << '(' << i << ',' << i << ") " << diagonal
444 << std::endl;
445 break;
446 }
447 }
448 }
449 for (size_type j = rows[i]; j < rows[i + 1]; ++j)
450 {
451 if (has_diagonal && i == cols[j])
452 continue;
453 if (across)
454 out << ' ' << i << ',' << cols[j] << ':' << val[j];
455 else
456 out << '(' << i << ',' << cols[j] << ") " << val[j] << std::endl;
457 }
458 }
459 if (across)
460 out << std::endl;
461 }
462
463
464
465 template <typename Number>
466 void
468 const unsigned int precision,
469 const bool scientific,
470 const unsigned int width_,
471 const char *zero_string,
472 const double denominator,
473 const char *separator) const
474 {
475 Assert(column_index_dev != nullptr, ExcNotInitialized());
476 Assert(val_dev != nullptr, ExcNotInitialized());
477 Assert(row_ptr_dev != nullptr, ExcNotInitialized());
478
479 std::vector<int> rows(n_rows + 1);
480 std::vector<int> cols(nnz);
481 std::vector<Number> val(nnz);
482 Utilities::CUDA::copy_to_host(row_ptr_dev.get(), rows);
483 Utilities::CUDA::copy_to_host(column_index_dev.get(), cols);
484 Utilities::CUDA::copy_to_host(val_dev.get(), val);
485
486 unsigned int width = width_;
487
488 std::ios::fmtflags old_flags = out.flags();
489 unsigned int old_precision = out.precision(precision);
490
491 if (scientific)
492 {
493 out.setf(std::ios::scientific, std::ios::floatfield);
494 if (!width)
495 width = precision + 7;
496 }
497 else
498 {
499 out.setf(std::ios::fixed, std::ios::floatfield);
500 if (!width)
501 width = precision + 2;
502 }
503
504 for (size_type i = 0; i < n_rows; ++i)
505 {
506 size_type j = rows[i];
507 for (size_type k = 0; k < n_cols; ++k)
508 {
509 if (k == cols[j])
510 {
511 out << std::setw(width) << val[j] * Number(denominator)
512 << separator;
513 ++j;
514 }
515 else
516 out << std::setw(width) << zero_string << separator;
517 }
518 out << std::endl;
519 };
520 AssertThrow(out.fail() == false, ExcIO());
521
522 // reset output format
523 out.precision(old_precision);
524 out.flags(old_flags);
525 }
526} // namespace CUDAWrappers
527
529
530#endif
531#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
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 char *separator=" ") const
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
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:131