42#ifndef TPETRA_BLOCKVIEW_HPP
43#define TPETRA_BLOCKVIEW_HPP
53#include "TpetraCore_config.h"
54#include "Kokkos_ArithTraits.hpp"
55#include "Kokkos_Complex.hpp"
72template<
class ViewType1,
74 const int rank1 = ViewType1::rank>
92 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
93 "AbsMax: ViewType1 and ViewType2 must have the same rank.");
94 typedef typename std::remove_reference<
decltype (
Y(0,0)) >::type
STY;
95 static_assert(! std::is_const<STY>::value,
96 "AbsMax: The type of each entry of Y must be nonconst.");
97 typedef typename std::decay<
decltype (
X(0,0)) >::type
STX;
98 static_assert( std::is_same<STX, STY>::value,
99 "AbsMax: The type of each entry of X and Y must be the same.");
100 typedef Kokkos::Details::ArithTraits<STY>
KAT;
102 const int numCols =
Y.extent (1);
104 for (
int j = 0;
j < numCols; ++
j) {
132 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
133 "AbsMax: ViewType1 and ViewType2 must have the same rank.");
135 typedef typename std::remove_reference<
decltype (
Y(0)) >::type
STY;
136 static_assert(! std::is_const<STY>::value,
137 "AbsMax: The type of each entry of Y must be nonconst.");
138 typedef typename std::remove_const<
typename std::remove_reference<
decltype (
X(0)) >::type>::type
STX;
139 static_assert( std::is_same<STX, STY>::value,
140 "AbsMax: The type of each entry of X and Y must be the same.");
141 typedef Kokkos::Details::ArithTraits<STY>
KAT;
164template<
class ViewType1,
class ViewType2, const
int rank = ViewType1::rank>
168 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
169 "absMax: ViewType1 and ViewType2 must have the same rank.");
181 const int rank = ViewType::rank>
200#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
234 using x_value_type =
typename std::decay<
decltype (*
x.data()) >::type;
237#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
249template<
class ViewType,
251 class IndexType = int,
252 const bool is_contiguous =
false,
253 const int rank = ViewType::rank>
299#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
312template<
class CoefficientType,
315 class IndexType = int,
316 const bool is_contiguous =
false,
317 const int rank = ViewType1::rank>
338 using Kokkos::Details::ArithTraits;
339 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
340 "AXPY: x and y must have the same rank.");
364 using Kokkos::Details::ArithTraits;
365 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
366 "AXPY: X and Y must have the same rank.");
390 using Kokkos::Details::ArithTraits;
391 static_assert (
static_cast<int> (ViewType1::rank) ==
static_cast<int> (ViewType2::rank),
392 "AXPY: x and y must have the same rank.");
395 using x_value_type =
typename std::decay<
decltype (*
x.data()) >::type;
396 using y_value_type =
typename std::decay<
decltype (*
y.data()) >::type;
400#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
413template<
class ViewType1,
415 class IndexType = int,
416 const bool is_contiguous =
false,
417 const int rank = ViewType1::rank>
469 using x_value_type =
typename std::decay<
decltype (*
x.data()) >::type;
470 using y_value_type =
typename std::decay<
decltype (*
y.data()) >::type;
474#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
482template<
class VecType1,
486 class IndexType = int,
487 bool is_contiguous =
false,
488 class BlkLayoutType =
typename BlkType::array_layout>
490 static KOKKOS_INLINE_FUNCTION
void
491 run (
const CoeffType& alpha,
496 static_assert (VecType1::rank == 1,
"GEMV: VecType1 must have rank 1.");
497 static_assert (BlkType::rank == 2,
"GEMV: BlkType must have rank 2.");
498 static_assert (VecType2::rank == 1,
"GEMV: VecType2 must have rank 1.");
500 const IndexType numRows =
static_cast<IndexType
> (A.extent (0));
501 const IndexType numCols =
static_cast<IndexType
> (A.extent (1));
504 for (IndexType i = 0; i < numRows; ++i)
505 for (IndexType j = 0; j < numCols; ++j)
506 y(i) += alpha * A(i,j) * x(j);
510template<
class VecType1,
515struct GEMV<VecType1,BlkType,VecType2,CoeffType,IndexType,true,Kokkos::LayoutLeft> {
516 static KOKKOS_INLINE_FUNCTION
void
517 run (
const CoeffType& alpha,
522 static_assert (VecType1::rank == 1,
"GEMV: VecType1 must have rank 1.");
523 static_assert (BlkType::rank == 2,
"GEMV: BlkType must have rank 2.");
524 static_assert (VecType2::rank == 1,
"GEMV: VecType2 must have rank 1.");
526 using A_value_type =
typename std::decay<
decltype (A(0,0)) >::type;
527 using x_value_type =
typename std::decay<
decltype (x(0)) >::type;
528 using y_value_type =
typename std::decay<
decltype (y(0)) >::type;
530 const IndexType numRows =
static_cast<IndexType
> (A.extent (0));
531 const IndexType numCols =
static_cast<IndexType
> (A.extent (1));
533 const A_value_type *__restrict__ A_ptr(A.data());
const IndexType as1(A.stride(1));
534 const x_value_type *__restrict__ x_ptr(x.data());
535 y_value_type *__restrict__ y_ptr(y.data());
537#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
540 for (IndexType j=0;j<numCols;++j) {
541 const x_value_type x_at_j = alpha*x_ptr[j];
542 const A_value_type *__restrict__ A_at_j = A_ptr + j*as1;
543#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
546 for (IndexType i=0;i<numRows;++i)
547 y_ptr[i] += A_at_j[i] * x_at_j;
552template<
class VecType1,
557struct GEMV<VecType1,BlkType,VecType2,CoeffType,IndexType,true,Kokkos::LayoutRight> {
558 static KOKKOS_INLINE_FUNCTION
void
559 run (
const CoeffType& alpha,
564 static_assert (VecType1::rank == 1,
"GEMV: VecType1 must have rank 1.");
565 static_assert (BlkType::rank == 2,
"GEMV: BlkType must have rank 2.");
566 static_assert (VecType2::rank == 1,
"GEMV: VecType2 must have rank 1.");
568 using A_value_type =
typename std::decay<
decltype (A(0,0)) >::type;
569 using x_value_type =
typename std::decay<
decltype (x(0)) >::type;
570 using y_value_type =
typename std::decay<
decltype (y(0)) >::type;
572 const IndexType numRows =
static_cast<IndexType
> (A.extent (0));
573 const IndexType numCols =
static_cast<IndexType
> (A.extent (1));
575 const A_value_type *__restrict__ A_ptr(A.data());
const IndexType as0(A.stride(0));
576 const x_value_type *__restrict__ x_ptr(x.data());
577 y_value_type *__restrict__ y_ptr(y.data());
579#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
582 for (IndexType i=0;i<numRows;++i) {
583 y_value_type y_at_i(0);
584 const auto A_at_i = A_ptr + i*as0;
585#if defined(KOKKOS_ENABLE_PRAGMA_UNROLL)
588 for (IndexType j=0;j<numCols;++j)
589 y_at_i += A_at_i[j] * x_ptr[j];
590 y_ptr[i] += alpha*y_at_i;
599template<
class ViewType,
600 class CoefficientType,
601 class IndexType = int,
602 const int rank = ViewType::rank>
603KOKKOS_INLINE_FUNCTION
void
606 using LayoutType =
typename ViewType::array_layout;
607 constexpr bool is_contiguous = (std::is_same<LayoutType,Kokkos::LayoutLeft>::value ||
608 std::is_same<LayoutType,Kokkos::LayoutRight>::value);
613template<
class ViewType,
615 class IndexType = int,
616 const int rank = ViewType::rank>
617KOKKOS_INLINE_FUNCTION
void
620 using LayoutType =
typename ViewType::array_layout;
621 constexpr bool is_contiguous = (std::is_same<LayoutType,Kokkos::LayoutLeft>::value ||
622 std::is_same<LayoutType,Kokkos::LayoutRight>::value);
631template<
class CoefficientType,
634 class IndexType = int,
635 const int rank = ViewType1::rank>
636KOKKOS_INLINE_FUNCTION
void
641 static_assert (
static_cast<int> (ViewType1::rank) ==
642 static_cast<int> (ViewType2::rank),
643 "AXPY: x and y must have the same rank.");
644 using LayoutType1 =
typename ViewType1::array_layout;
645 using LayoutType2 =
typename ViewType2::array_layout;
646 constexpr bool is_layout_same = std::is_same<LayoutType1,LayoutType2>::value;
647 constexpr bool is_x_contiguous = (std::is_same<LayoutType1,Kokkos::LayoutLeft>::value ||
648 std::is_same<LayoutType1,Kokkos::LayoutRight>::value);
661template<
class ViewType1,
663 class IndexType = int,
664 const int rank = ViewType1::rank>
665KOKKOS_INLINE_FUNCTION
void
668 static_assert (
static_cast<int> (ViewType1::rank) ==
669 static_cast<int> (ViewType2::rank),
670 "COPY: x and y must have the same rank.");
671 using LayoutType1 =
typename ViewType1::array_layout;
672 using LayoutType2 =
typename ViewType2::array_layout;
673 constexpr bool is_layout_same = std::is_same<LayoutType1,LayoutType2>::value;
674 constexpr bool is_x_contiguous = (std::is_same<LayoutType1,Kokkos::LayoutLeft>::value ||
675 std::is_same<LayoutType1,Kokkos::LayoutRight>::value);
691template<
class VecType1,
695 class IndexType =
int>
696KOKKOS_INLINE_FUNCTION
void
702 constexpr bool is_A_contiguous = (std::is_same<typename BlkType::array_layout, Kokkos::LayoutLeft>::value ||
703 std::is_same<typename BlkType::array_layout, Kokkos::LayoutRight>::value);
704 constexpr bool is_x_contiguous = (std::is_same<typename VecType1::array_layout, Kokkos::LayoutLeft>::value ||
705 std::is_same<typename VecType1::array_layout, Kokkos::LayoutRight>::value);
706 constexpr bool is_y_contiguous = (std::is_same<typename VecType2::array_layout, Kokkos::LayoutLeft>::value ||
707 std::is_same<typename VecType2::array_layout, Kokkos::LayoutRight>::value);
709 Impl::GEMV<VecType1, BlkType, VecType2, CoeffType, IndexType, is_contiguous>::run (
alpha,
A,
x,
y);
719template<
class ViewType1,
722 class CoefficientType,
723 class IndexType =
int>
724KOKKOS_INLINE_FUNCTION
void
734 static_assert (ViewType1::rank == 2,
"GEMM: A must have rank 2 (be a matrix).");
735 static_assert (ViewType2::rank == 2,
"GEMM: B must have rank 2 (be a matrix).");
736 static_assert (ViewType3::rank == 2,
"GEMM: C must have rank 2 (be a matrix).");
738 typedef typename std::remove_reference<
decltype (
A(0,0))>::type
Scalar;
739 typedef Kokkos::Details::ArithTraits<Scalar> STS;
861template<
class LittleBlockType,
862 class LittleVectorType>
863KOKKOS_INLINE_FUNCTION
void
867 typedef typename std::decay<
decltype (
ipiv(0)) >::type
IndexType;
868 static_assert (std::is_integral<IndexType>::value,
869 "GETF2: The type of each entry of ipiv must be an integer type.");
870 typedef typename std::remove_reference<
decltype (
A(0,0))>::type
Scalar;
871 static_assert (! std::is_const<Scalar>::value,
872 "GETF2: A must not be a const View (or LittleBlock).");
873 static_assert (! std::is_const<std::remove_reference<
decltype (
ipiv(0))>>::value,
874 "GETF2: ipiv must not be a const View (or LittleBlock).");
875 static_assert (LittleBlockType::rank == 2,
"GETF2: A must have rank 2 (be a matrix).");
876 typedef Kokkos::Details::ArithTraits<Scalar> STS;
899 if(STS::abs(
A(
i,
j)) > STS::abs(
A(
jp,
j))) {
942template<
class LittleBlockType,
943 class LittleIntVectorType,
944 class LittleScalarVectorType,
945 const int rank = LittleScalarVectorType::rank>
948 run (
const char mode[],
961 run (
const char mode[],
968 typedef typename std::decay<
decltype (
ipiv(0))>::type
IndexType;
971 static_assert (std::is_integral<IndexType>::value &&
972 std::is_signed<IndexType>::value,
973 "GETRS: The type of each entry of ipiv must be a signed integer.");
974 typedef typename std::decay<
decltype (
A(0,0))>::type
Scalar;
975 static_assert (! std::is_const<std::remove_reference<
decltype (
B(0))>>::value,
976 "GETRS: B must not be a const View (or LittleBlock).");
977 static_assert (LittleBlockType::rank == 2,
"GETRS: A must have rank 2 (be a matrix).");
978 static_assert (LittleIntVectorType::rank == 1,
"GETRS: ipiv must have rank 1.");
979 static_assert (LittleScalarVectorType::rank == 1,
"GETRS: For this specialization, B must have rank 1.");
981 typedef Kokkos::Details::ArithTraits<Scalar> STS;
1002 if(
mode[0] ==
'n' ||
mode[0] ==
'N') {
1034 else if(
mode[0] ==
't' ||
mode[0] ==
'T') {
1039 else if (
mode[0] ==
'c' ||
mode[0] ==
'C') {
1057 run (
const char mode[],
1064 typedef typename std::decay<
decltype (
ipiv(0)) >::type
IndexType;
1065 static_assert (std::is_integral<IndexType>::value,
1066 "GETRS: The type of each entry of ipiv must be an integer type.");
1067 static_assert (! std::is_const<std::remove_reference<
decltype (
B(0)) > >::value,
1068 "GETRS: B must not be a const View (or LittleBlock).");
1069 static_assert (LittleBlockType::rank == 2,
"GETRS: A must have rank 2 (be a matrix).");
1070 static_assert (LittleIntVectorType::rank == 1,
"GETRS: ipiv must have rank 1.");
1071 static_assert (LittleScalarVectorType::rank == 2,
"GETRS: For this specialization, B must have rank 2.");
1080 auto B_cur = Kokkos::subview (
B, Kokkos::ALL (),
rhs);
1123template<
class LittleBlockType,
1124 class LittleIntVectorType,
1125 class LittleScalarVectorType>
1126KOKKOS_INLINE_FUNCTION
void
1133 typedef typename std::decay<
decltype (
ipiv(0))>::type
IndexType;
1136 static_assert (std::is_integral<IndexType>::value &&
1137 std::is_signed<IndexType>::value,
1138 "GETRI: The type of each entry of ipiv must be a signed integer.");
1139 typedef typename std::remove_reference<
decltype (
A(0,0))>::type
Scalar;
1140 static_assert (! std::is_const<std::remove_reference<
decltype (
A(0,0))>>::value,
1141 "GETRI: A must not be a const View (or LittleBlock).");
1142 static_assert (! std::is_const<std::remove_reference<
decltype (
work(0))>>::value,
1143 "GETRI: work must not be a const View (or LittleBlock).");
1144 static_assert (LittleBlockType::rank == 2,
"GETRI: A must have rank 2 (be a matrix).");
1145 typedef Kokkos::Details::ArithTraits<Scalar> STS;
1228template<
class LittleBlockType,
1229 class LittleVectorTypeX,
1230 class LittleVectorTypeY,
1231 class CoefficientType,
1232 class IndexType =
int>
1233KOKKOS_INLINE_FUNCTION
void
1234GEMV (
const char trans,
1235 const CoefficientType& alpha,
1236 const LittleBlockType& A,
1237 const LittleVectorTypeX& x,
1238 const CoefficientType& beta,
1239 const LittleVectorTypeY& y)
1245 typedef typename std::remove_reference<
decltype (y(0)) >::type y_value_type;
1246 const IndexType numRows =
static_cast<IndexType
> (A.extent (0));
1247 const IndexType numCols =
static_cast<IndexType
> (A.extent (1));
1251 for (IndexType i = 0; i < numRows; ++i) {
1256 for (IndexType i = 0; i < numRows; ++i) {
1257 y_value_type y_i = 0.0;
1258 for (IndexType j = 0; j < numCols; ++j) {
1259 y_i += A(i,j) * x(j);
1268 for (IndexType i = 0; i < numRows; ++i) {
1273 for (IndexType i = 0; i < numRows; ++i) {
1279 for (IndexType i = 0; i < numRows; ++i) {
1280 y_value_type y_i = beta * y(i);
1281 for (IndexType j = 0; j < numCols; ++j) {
1282 y_i += alpha * A(i,j) * x(j);
Struct that holds views of the contents of a CrsMatrix.
KOKKOS_INLINE_FUNCTION void absMax(const ViewType2 &Y, const ViewType1 &X)
Implementation of Tpetra's ABSMAX CombineMode for the small dense blocks in BlockCrsMatrix,...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
KOKKOS_INLINE_FUNCTION void GEMV(const CoeffType &alpha, const BlkType &A, const VecType1 &x, const VecType2 &y)
y := y + alpha * A * x (dense matrix-vector multiply)
KOKKOS_INLINE_FUNCTION void GETF2(const LittleBlockType &A, const LittleVectorType &ipiv, int &info)
Computes A = P*L*U.
KOKKOS_INLINE_FUNCTION void GETRI(const LittleBlockType &A, const LittleIntVectorType &ipiv, const LittleScalarVectorType &work, int &info)
Compute inverse of A, using result of GETRF or GETF2.
KOKKOS_INLINE_FUNCTION void AXPY(const CoefficientType &alpha, const ViewType1 &x, const ViewType2 &y)
y := y + alpha * x (dense vector or matrix update)
KOKKOS_INLINE_FUNCTION void FILL(const ViewType &x, const InputType &val)
Set every entry of x to val.
KOKKOS_INLINE_FUNCTION void SCAL(const CoefficientType &alpha, const ViewType &x)
x := alpha*x, where x is either rank 1 (a vector) or rank 2 (a matrix).
KOKKOS_INLINE_FUNCTION void COPY(const ViewType1 &x, const ViewType2 &y)
Deep copy x into y, where x and y are either rank 1 (vectors) or rank 2 (matrices) with the same dime...
KOKKOS_INLINE_FUNCTION void GETRS(const char mode[], const LittleBlockType &A, const LittleIntVectorType &ipiv, const LittleScalarVectorType &B, int &info)
Solve the linear system(s) AX=B, using the result of GETRF or GETF2.
KOKKOS_INLINE_FUNCTION void GEMM(const char transA[], const char transB[], const CoefficientType &alpha, const ViewType1 &A, const ViewType2 &B, const CoefficientType &beta, const ViewType3 &C)
Small dense matrix-matrix multiply: C := alpha*A*B + beta*C
@ ZERO
Replace old values with zero.
static KOKKOS_INLINE_FUNCTION void run(const CoefficientType &alpha, const ViewType1 &x, const ViewType2 &y)
y := y + alpha*x (rank-1 x and y, i.e., vectors)
static KOKKOS_INLINE_FUNCTION void run(const CoefficientType &alpha, const ViewType1 &X, const ViewType2 &Y)
Y := Y + alpha*X (rank-2 X and Y, i.e., matrices)
Implementation of Tpetra::AXPY function.
static KOKKOS_INLINE_FUNCTION void run(const ViewType2 &Y, const ViewType1 &X)
(*this)(i) := max(abs((*this)(i)), abs(X(i))) for all i.
static KOKKOS_INLINE_FUNCTION void run(const ViewType2 &Y, const ViewType1 &X)
(*this)(i,j) := max(abs((*this)(i,j)), abs(X(i,j))) for all (i,j).
Implementation of Tpetra's ABSMAX CombineMode for the small dense blocks in BlockCrsMatrix,...
static KOKKOS_INLINE_FUNCTION void run(const ViewType1 &x, const ViewType2 &y)
y := x (rank-1 x and y, i.e., vectors)
static KOKKOS_INLINE_FUNCTION void run(const ViewType1 &X, const ViewType2 &Y)
Y := X (rank-2 X and Y, i.e., matrices)
Implementation of Tpetra::COPY function.
Implementation of Tpetra::FILL function.
Computes the solution to Ax=b.
static KOKKOS_INLINE_FUNCTION void run(const CoefficientType &alpha, const ViewType &x)
x := alpha*x (rank-1 x, i.e., a vector)
static KOKKOS_INLINE_FUNCTION void run(const CoefficientType &alpha, const ViewType &A)
A := alpha*A (rank-2 A, i.e., a matrix)
Implementation of Tpetra::SCAL function.