Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
57template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
58class MatVec {
59 const ViewTypeA A;
60 const ViewTypeB b;
61 const ViewTypeC c;
62 const size_t m, n;
63
64public:
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
74 KOKKOS_INLINE_FUNCTION
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
88template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
89void 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.
96template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
97class MatVecDeriv {
98 const ViewTypeA A;
99 const ViewTypeB b;
100 const ViewTypeC c;
101 const size_t m, n, p;
102
103public:
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
113 KOKKOS_INLINE_FUNCTION
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
134template <typename ViewTypeA, typename ViewTypeB, typename ViewTypeC>
135void 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
143int 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}
expr expr expr fastAccessDx(i)) FAD_UNARYOP_MACRO(exp
expr val()
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
Sacado::Fad::DFad< double > FadType
void run_mat_vec(const ViewTypeA &A, const ViewTypeB &b, const ViewTypeC &c)
char c_
const char * p
void run_mat_vec_deriv(const ViewTypeA &A, const ViewTypeB &b, const ViewTypeC &c)
Definition: mat_vec.cpp:98
const double tol