Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
blas_example.cpp
Go to the documentation of this file.
1 // $Id$
2 // $Source$
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Sacado Package
7 // Copyright (2006) Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // This library is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as
14 // published by the Free Software Foundation; either version 2.1 of the
15 // License, or (at your option) any later version.
16 //
17 // This library is distributed in the hope that it will be useful, but
18 // WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 // Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License along with this library; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
25 // USA
26 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
27 // (etphipp@sandia.gov).
28 //
29 // ***********************************************************************
30 // @HEADER
31 
32 // blas_example
33 //
34 // usage:
35 // blas_example
36 //
37 // output:
38 // prints the results of differentiating a BLAS routine with forward
39 // mode AD using the Sacado::Fad::DFad class (uses dynamic memory
40 // allocation for number of derivative components).
41 
42 #include <iostream>
43 #include <iomanip>
44 
45 #include "Sacado_No_Kokkos.hpp"
46 #include "Teuchos_BLAS.hpp"
47 #include "Sacado_Fad_BLAS.hpp"
48 
50 
51 int main(int argc, char **argv)
52 {
53  const unsigned int n = 5;
54  std::vector<double> a(n*n), b(n), c(n);
55  std::vector<FadType> A(n*n), B(n), C(n);
56  for (unsigned int i=0; i<n; i++) {
57  for (unsigned int j=0; j<n; j++)
60  c[i] = 0.0;
61 
62  for (unsigned int j=0; j<n; j++)
63  A[i+j*n] = FadType(a[i+j*n]);
64  B[i] = FadType(n, i, b[i]);
65  C[i] = FadType(c[i]);
66  }
67 
69  blas.GEMV(Teuchos::NO_TRANS, n, n, 1.0, &a[0], n, &b[0], 1, 0.0, &c[0], 1);
70 
71  // Teuchos::BLAS<int,FadType> fad_blas;
72  // fad_blas.GEMV(Teuchos::NO_TRANS, n, n, 1.0, &A[0], n, &B[0], 1, 0.0, &C[0], 1);
73 
74  Teuchos::BLAS<int,FadType> sacado_fad_blas(false,false,3*n*n+2*n);
75  sacado_fad_blas.GEMV(Teuchos::NO_TRANS, n, n, 1.0, &A[0], n, &B[0], 1, 0.0, &C[0], 1);
76 
77  // Print the results
78  int p = 4;
79  int w = p+7;
80  std::cout.setf(std::ios::scientific);
81  std::cout.precision(p);
82 
83  std::cout << "BLAS GEMV calculation:" << std::endl;
84  std::cout << "a = " << std::endl;
85  for (unsigned int i=0; i<n; i++) {
86  for (unsigned int j=0; j<n; j++)
87  std::cout << " " << std::setw(w) << a[i+j*n];
88  std::cout << std::endl;
89  }
90  std::cout << "b = " << std::endl;
91  for (unsigned int i=0; i<n; i++) {
92  std::cout << " " << std::setw(w) << b[i];
93  }
94  std::cout << std::endl;
95  std::cout << "c = " << std::endl;
96  for (unsigned int i=0; i<n; i++) {
97  std::cout << " " << std::setw(w) << c[i];
98  }
99  std::cout << std::endl << std::endl;
100 
101  std::cout << "FAD BLAS GEMV calculation:" << std::endl;
102  std::cout << "A.val() (should = a) = " << std::endl;
103  for (unsigned int i=0; i<n; i++) {
104  for (unsigned int j=0; j<n; j++)
105  std::cout << " " << std::setw(w) << A[i+j*n].val();
106  std::cout << std::endl;
107  }
108  std::cout << "B.val() (should = b) = " << std::endl;
109  for (unsigned int i=0; i<n; i++) {
110  std::cout << " " << std::setw(w) << B[i].val();
111  }
112  std::cout << std::endl;
113  std::cout << "C.val() (should = c) = " << std::endl;
114  for (unsigned int i=0; i<n; i++) {
115  std::cout << " " << std::setw(w) << C[i].val();
116  }
117  std::cout << std::endl;
118  std::cout << "C.dx() ( = dc/db, should = a) = " << std::endl;
119  for (unsigned int i=0; i<n; i++) {
120  for (unsigned int j=0; j<n; j++)
121  std::cout << " " << std::setw(w) << C[i].dx(j);
122  std::cout << std::endl;
123  }
124 
125  double tol = 1.0e-14;
126  bool failed = false;
127  for (unsigned int i=0; i<n; i++) {
128  if (std::fabs(C[i].val() - c[i]) > tol)
129  failed = true;
130  for (unsigned int j=0; j<n; j++) {
131  if (std::fabs(C[i].dx(j) - a[i+j*n]) > tol)
132  failed = true;
133  }
134  }
135  if (!failed) {
136  std::cout << "\nExample passed!" << std::endl;
137  return 0;
138  }
139  else {
140  std::cout <<"\nSomething is wrong, example failed!" << std::endl;
141  return 1;
142  }
143 }
expr expr dx(i)
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const
Sacado::Fad::DFad< double > FadType
expr val()
#define C(x)
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
const double tol
void GEMV(ETransp trans, const int &m, const int &n, const double &alpha, const double *A, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy) const
fabs(expr.val())
int n