Sacado Package Browser (Single Doxygen Collection)  Version of the Day
view_example.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include "Sacado.hpp"
43 
44 //
45 // A simple example showing how to use Sacado with Kokkos
46 //
47 // The code would be simpler using lambda's instead of functors, but that
48 // can be problematic for Cuda and for versions of g++ that claim to, but
49 // don't really support C++11
50 //
51 
52 // This code relies on the Kokkos view specializations being enabled, which
53 // is the default, but can be disabled by the user
54 #if defined(HAVE_SACADO_VIEW_SPEC) && !defined(SACADO_DISABLE_FAD_VIEW_SPEC)
55 
56 // Functor to compute matrix-vector product c = A*b using Kokkos
57 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
58 class MatVec {
59  const ViewTypeA A;
60  const ViewTypeB b;
61  const ViewTypeC c;
62  const size_t m, n;
63 
64 public:
65 
66  MatVec(const ViewTypeA& A_, const ViewTypeB& b_, const ViewTypeC& c_) :
67  A(A_), b(b_), c(c_),
68  m(A.extent(0)), n(A.extent(1))
69  {
70  typedef typename ViewTypeC::execution_space execution_space;
71  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>(0,m), *this);
72  }
73 
75  void operator() (const size_t i) const {
76  typedef typename ViewTypeC::value_type scalar_type;
77 
78  scalar_type t = 0.0;
79  for (size_t j=0; j<n; ++j)
80  t += A(i,j)*b(j);
81  c(i) = t;
82  }
83 };
84 
85 // Function to run mat-vec functor above
86 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
87 void run_mat_vec(const ViewTypeA& A, const ViewTypeB& b,const ViewTypeC& c)
88 {
89  MatVec<ViewTypeA,ViewTypeB,ViewTypeC> f(A,b,c);
90 }
91 
92 // Functor to compute the derivative of the matrix-vector product
93 // c = A*b using Kokkos.
94 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
95 class MatVecDeriv {
96  const ViewTypeA A;
97  const ViewTypeB b;
98  const ViewTypeC c;
99  const size_t m, n, p;
100 
101 public:
102 
103  MatVecDeriv(const ViewTypeA& A_, const ViewTypeB& b_, const ViewTypeC& c_) :
104  A(A_), b(b_), c(c_),
105  m(A.extent(0)), n(A.extent(1)), p(A.extent(2)-1)
106  {
107  typedef typename ViewTypeC::execution_space execution_space;
108  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>(0,m), *this);
109  }
110 
112  void operator() (const size_t i) const {
113  typedef typename ViewTypeC::value_type scalar_type;
114 
115  // Derivative portion
116  for (size_t k=0; k<p; ++k) {
117  scalar_type t = 0.0;
118  for (size_t j=0; j<n; ++j)
119  t += A(i,j,k)*b(j,p) + A(i,j,p)*b(j,k);
120  c(i,k) = t;
121  }
122 
123  // Value portion
124  scalar_type t = 0.0;
125  for (size_t j=0; j<n; ++j)
126  t += A(i,j,p)*b(j,p);
127  c(i,p) = t;
128  }
129 };
130 
131 // Function to run mat-vec derivative functor above
132 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
133 void run_mat_vec_deriv(const ViewTypeA& A, const ViewTypeB& b,
134  const ViewTypeC& c)
135 {
136  MatVecDeriv<ViewTypeA,ViewTypeB,ViewTypeC> f(A,b,c);
137 }
138 
139 #endif
140 
141 int main(int argc, char* argv[]) {
142  int ret = 0;
143 
144 #if defined(HAVE_SACADO_VIEW_SPEC) && !defined(SACADO_DISABLE_FAD_VIEW_SPEC)
145 
146  Kokkos::initialize(argc, argv);
147  {
148 
149  const size_t m = 10; // Number of rows
150  const size_t n = 2; // Number of columns
151  const size_t p = 1; // Derivative dimension
152 
153  // Allocate Kokkos view's for matrix, input vector, and output vector
154  // Derivative dimension must be specified as the last constructor argument
156  Kokkos::View<FadType**> A("A",m,n,p+1);
157  Kokkos::View<FadType*> b("b",n,p+1);
158  Kokkos::View<FadType*> c("c",m,p+1);
159 
160  // Initialize values
161  Kokkos::deep_copy( A, FadType(p, 0, 2.0) );
162  Kokkos::deep_copy( b, FadType(p, 0, 3.0) );
163  Kokkos::deep_copy( c, 0.0 );
164 
165  // Run mat-vec
166  run_mat_vec(A,b,c);
167 
168  // Print result
169  std::cout << "\nc = A*b: Differentiated using Sacado:" << std::endl;
170  auto h_c = Kokkos::create_mirror_view(c);
171  Kokkos::deep_copy(h_c, c);
172  for (size_t i=0; i<m; ++i)
173  std::cout << "\tc(" << i << ") = " << h_c(i) << std::endl;
174 
175  // Now compute derivative analytically. Any Sacado view can be flattened
176  // into a standard view of one higher rank, with the extra dimension equal
177  // to p+1
178  Kokkos::View<FadType*> c2("c",m,p+1);
179  Kokkos::View<double***> A_flat = A;
180  Kokkos::View<double**> b_flat = b;
181  Kokkos::View<double**> c_flat = c2;
182  run_mat_vec_deriv(A_flat, b_flat, c_flat);
183 
184  // Print result
185  std::cout << "\nc = A*b: Differentiated analytically:" << std::endl;
186  auto h_c2 = Kokkos::create_mirror_view(c2);
187  Kokkos::deep_copy(h_c2, c2);
188  for (size_t i=0; i<m; ++i)
189  std::cout << "\tc(" << i << ") = " << h_c2(i) << std::endl;
190 
191  // Compute the error
192  double err = 0.0;
193  for (size_t i=0; i<m; ++i) {
194  for (size_t k=0; k<p; ++k)
195  err += std::abs(h_c(i).fastAccessDx(k)-h_c2(i).fastAccessDx(k));
196  err += std::abs(h_c(i).val()-h_c2(i).val());
197  }
198 
199  double tol = 1.0e-14;
200  if (err < tol) {
201  ret = 0;
202  std::cout << "\nExample passed!" << std::endl;
203  }
204  else {
205  ret = 1;
206  std::cout <<"\nSomething is wrong, example failed!" << std::endl;
207  }
208 
209  }
210  Kokkos::finalize();
211 
212 #endif
213 
214  return ret;
215 }
abs(expr.val())
void run_mat_vec(const ViewTypeA &A, const ViewTypeB &b, const ViewTypeC &c)
Sacado::Fad::DFad< double > FadType
expr val()
#define KOKKOS_INLINE_FUNCTION
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
#define A
Definition: Sacado_rad.hpp:572
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
int main(int argc, char *argv[])
const double tol
void run_mat_vec_deriv(const ViewTypeA &A, const ViewTypeB &b, const ViewTypeC &c)