Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  scalar_type bb = b(j); // fix for Intel 17.0.1 with optimization
81  t += A(i,j)*bb;
82  }
83  c(i) = t;
84  }
85 };
86 
87 // Function to run mat-vec functor above
88 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
89 void run_mat_vec(const ViewTypeA& A, const ViewTypeB& b,const ViewTypeC& c)
90 {
91  MatVec<ViewTypeA,ViewTypeB,ViewTypeC> f(A,b,c);
92 }
93 
94 // Functor to compute the derivative of the matrix-vector product
95 // c = A*b using Kokkos.
96 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
97 class MatVecDeriv {
98  const ViewTypeA A;
99  const ViewTypeB b;
100  const ViewTypeC c;
101  const size_t m, n, p;
102 
103 public:
104 
105  MatVecDeriv(const ViewTypeA& A_, const ViewTypeB& b_, const ViewTypeC& c_) :
106  A(A_), b(b_), c(c_),
107  m(A.extent(0)), n(A.extent(1)), p(A.extent(2)-1)
108  {
109  typedef typename ViewTypeC::execution_space execution_space;
110  Kokkos::parallel_for(Kokkos::RangePolicy<execution_space>(0,m), *this);
111  }
112 
114  void operator() (const size_t i) const {
115  typedef typename ViewTypeC::value_type scalar_type;
116 
117  // Derivative portion
118  for (size_t k=0; k<p; ++k) {
119  scalar_type t = 0.0;
120  for (size_t j=0; j<n; ++j)
121  t += A(i,j,k)*b(j,p) + A(i,j,p)*b(j,k);
122  c(i,k) = t;
123  }
124 
125  // Value portion
126  scalar_type t = 0.0;
127  for (size_t j=0; j<n; ++j)
128  t += A(i,j,p)*b(j,p);
129  c(i,p) = t;
130  }
131 };
132 
133 // Function to run mat-vec derivative functor above
134 template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
135 void run_mat_vec_deriv(const ViewTypeA& A, const ViewTypeB& b,
136  const ViewTypeC& c)
137 {
138  MatVecDeriv<ViewTypeA,ViewTypeB,ViewTypeC> f(A,b,c);
139 }
140 
141 #endif
142 
143 int main(int argc, char* argv[]) {
144  int ret = 0;
145 
146 #if defined(HAVE_SACADO_VIEW_SPEC) && !defined(SACADO_DISABLE_FAD_VIEW_SPEC)
147 
148  Kokkos::initialize(argc, argv);
149  {
150 
151  const size_t m = 10; // Number of rows
152  const size_t n = 2; // Number of columns
153  const size_t p = 1; // Derivative dimension
154 
155  // Allocate Kokkos view's for matrix, input vector, and output vector
156  // Derivative dimension must be specified as the last constructor argument
158  Kokkos::View<FadType**> A("A",m,n,p+1);
159  Kokkos::View<FadType*> b("b",n,p+1);
160  Kokkos::View<FadType*> c("c",m,p+1);
161 
162  // Initialize values
163  Kokkos::deep_copy( A, FadType(p, 0, 2.0) );
164  Kokkos::deep_copy( b, FadType(p, 0, 3.0) );
165  Kokkos::deep_copy( c, 0.0 );
166 
167  // Run mat-vec
168  run_mat_vec(A,b,c);
169 
170  // Print result
171  std::cout << "\nc = A*b: Differentiated using Sacado:" << std::endl;
172  auto h_c = Kokkos::create_mirror_view(c);
173  Kokkos::deep_copy(h_c, c);
174  for (size_t i=0; i<m; ++i)
175  std::cout << "\tc(" << i << ") = " << h_c(i) << std::endl;
176 
177  // Now compute derivative analytically. Any Sacado view can be flattened
178  // into a standard view of one higher rank, with the extra dimension equal
179  // to p+1
180  Kokkos::View<FadType*> c2("c",m,p+1);
181  Kokkos::View<double***> A_flat = A;
182  Kokkos::View<double**> b_flat = b;
183  Kokkos::View<double**> c_flat = c2;
184  run_mat_vec_deriv(A_flat, b_flat, c_flat);
185 
186  // Print result
187  std::cout << "\nc = A*b: Differentiated analytically:" << std::endl;
188  auto h_c2 = Kokkos::create_mirror_view(c2);
189  Kokkos::deep_copy(h_c2, c2);
190  for (size_t i=0; i<m; ++i)
191  std::cout << "\tc(" << i << ") = " << h_c2(i) << std::endl;
192 
193  // Compute the error
194  double err = 0.0;
195  for (size_t i=0; i<m; ++i) {
196  for (size_t k=0; k<p; ++k)
197  err += std::abs(h_c(i).fastAccessDx(k)-h_c2(i).fastAccessDx(k));
198  err += std::abs(h_c(i).val()-h_c2(i).val());
199  }
200 
201  double tol = 1.0e-14;
202  if (err < tol) {
203  ret = 0;
204  std::cout << "\nExample passed!" << std::endl;
205  }
206  else {
207  ret = 1;
208  std::cout <<"\nSomething is wrong, example failed!" << std::endl;
209  }
210 
211  }
212  Kokkos::finalize();
213 
214 #endif
215 
216  return ret;
217 }
void f()
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
int main()
Definition: ad_example.cpp:191
void run_mat_vec_deriv(const ViewTypeA &A, const ViewTypeB &b, const ViewTypeC &c)
Definition: mat_vec.cpp:98
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
const double tol
int n