Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
51int 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++)
58 a[i+j*n] = Teuchos::ScalarTraits<double>::random();
59 b[i] = Teuchos::ScalarTraits<double>::random();
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
68 Teuchos::BLAS<int,double> blas;
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)
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
#define C(x)
int main()
Sacado::Fad::DFad< double > FadType
const char * p
const double tol