Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_BLAS_UQ_PCE.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Teuchos: Common Tools Package
5// Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef _TEUCHOS_BLAS_UQ_PCE_HPP_
43#define _TEUCHOS_BLAS_UQ_PCE_HPP_
44
45#include "Teuchos_BLAS.hpp"
46#include "Sacado_UQ_PCE.hpp"
47
48// Specialize some things used in the default BLAS implementation that
49// don't seem correct for UQ::PCE scalar type
50namespace Teuchos {
51
52 namespace details {
53
54 template<typename Storage>
55 class GivensRotator<Sacado::UQ::PCE<Storage>, false> {
56 public:
59
60 void
62 ScalarType* db,
63 ScalarType* c,
64 ScalarType* s) const {
65
66 typedef ScalarTraits<ScalarType> STS;
67
68 // This is a straightforward translation into C++ of the
69 // reference BLAS' implementation of DROTG. You can get
70 // the Fortran 77 source code of DROTG here:
71 //
72 // http://www.netlib.org/blas/drotg.f
73 //
74 // I used the following rules to translate Fortran types and
75 // intrinsic functions into C++:
76 //
77 // DOUBLE PRECISION -> ScalarType
78 // DABS -> STS::magnitude
79 // DSQRT -> STM::squareroot
80 // DSIGN -> SIGN (see below)
81 //
82 // DSIGN(x,y) (the old DOUBLE PRECISION type-specific form of
83 // the Fortran type-generic SIGN intrinsic) required special
84 // translation, which we did in a separate utility function in
85 // the specializaton of GivensRotator for real arithmetic.
86 // (ROTG for complex arithmetic doesn't require this function.)
87 // C99 provides a copysign() math library function, but we are
88 // not able to rely on the existence of C99 functions here.
89 ScalarType r, roe, scale, z;
90
91 roe = *db;
92 if (STS::magnitude (*da) > STS::magnitude (*db)) {
93 roe = *da;
94 }
95 scale = STS::magnitude (*da) + STS::magnitude (*db);
96 if (scale == STS::zero()) {
97 *c = STS::one();
98 *s = STS::zero();
99 r = STS::zero();
100 z = STS::zero();
101 } else {
102 // I introduced temporaries into the translated BLAS code in
103 // order to make the expression easier to read and also save
104 // a few floating-point operations.
105 const ScalarType da_scaled = *da / scale;
106 const ScalarType db_scaled = *db / scale;
107 r = scale * STS::squareroot (da_scaled*da_scaled + db_scaled*db_scaled);
108 r = SIGN (STS::one(), roe) * r;
109 *c = *da / r;
110 *s = *db / r;
111 z = STS::one();
112 if (STS::magnitude (*da) > STS::magnitude (*db)) {
113 z = *s;
114 }
115 if (STS::magnitude (*db) >= STS::magnitude (*da) && *c != STS::zero()) {
116 z = STS::one() / *c;
117 }
118 }
119
120 *da = r;
121 *db = z;
122 }
123
124 private:
125
128 typedef typename ScalarType::value_type value_type;
129 typedef typename ScalarType::ordinal_type ordinal_type;
130
131 GivensRotator<value_type> value_rotator;
132 const ordinal_type sz = x.size() > y.size() ? x.size() : y.size();
133 ScalarType z(sz, 0.0);
134 for (ordinal_type i=0; i<sz; ++i)
135 z.fastAccessCoeff(i) = value_rotator.SIGN(x.coeff(i), y.coeff(i));
136 return z;
137 }
138 };
139
140 } // namespace details
141
142} // namespace Teuchos
143
144#endif // _TEUCHOS_BLAS_UQ_PCE_HPP_
void ROTG(ScalarType *da, ScalarType *db, ScalarType *c, ScalarType *s) const
ScalarType SIGN(ScalarType x, ScalarType y) const
Return ABS(x) if y > 0 or y is +0, else -ABS(x) (if y is -0 or < 0).