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_POLYNOMIAL_HPP
43 #define TEUCHOS_POLYNOMIAL_HPP
44
45 #include "Teuchos_PolynomialDecl.hpp"
46 #include "Teuchos_ScalarTraits.hpp"
47
48 template <typename CoeffT>
Polynomial(unsigned int deg,const CoeffT & cloneCoeff,unsigned int reserve)49 Teuchos::Polynomial<CoeffT>::Polynomial(unsigned int deg,
50 const CoeffT& cloneCoeff,
51 unsigned int reserve) :
52 d(deg)
53 {
54 if (reserve > d)
55 sz = reserve+1;
56 else
57 sz = d+1;
58
59 coeff.resize(sz);
60 for (unsigned int i=0; i<sz; i++)
61 coeff[i] = PolynomialTraits<CoeffT>::clone(cloneCoeff);
62 }
63
64 template <typename CoeffT>
Polynomial(unsigned int deg,unsigned int reserve)65 Teuchos::Polynomial<CoeffT>::Polynomial(unsigned int deg,
66 unsigned int reserve) :
67 d(deg)
68 {
69 if (reserve > d)
70 sz = reserve+1;
71 else
72 sz = d+1;
73
74 coeff.resize(sz);
75 }
76
77 template <typename CoeffT>
~Polynomial()78 Teuchos::Polynomial<CoeffT>::~Polynomial()
79 {
80 }
81
82 template <typename CoeffT>
83 void
setDegree(unsigned int deg)84 Teuchos::Polynomial<CoeffT>::setDegree(unsigned int deg)
85 {
86 d = deg;
87 if (d+1 > sz) {
88 coeff.resize(d+1);
89 if (coeff[0] != Teuchos::null) {
90 for (unsigned int i=sz; i<d+1; i++)
91 coeff[i] = PolynomialTraits<CoeffT>::clone(*coeff[0]);
92 }
93 sz = d+1;
94 }
95 }
96
97 template <typename CoeffT>
98 Teuchos::RCP<CoeffT>
getCoefficient(unsigned int i)99 Teuchos::Polynomial<CoeffT>::getCoefficient(unsigned int i)
100 {
101 #ifdef TEUCHOS_DEBUG
102 TEUCHOS_TEST_FOR_EXCEPTION(i > d,
103 std::out_of_range,
104 "Polynomial<CoeffT>::getCoefficient(i): " <<
105 "Error, coefficient i = " << i <<
106 " is not in range, degree = " << d << "." );
107 #endif
108 return coeff[i];
109 }
110
111 template <typename CoeffT>
112 Teuchos::RCP<const CoeffT>
getCoefficient(unsigned int i) const113 Teuchos::Polynomial<CoeffT>::getCoefficient(unsigned int i) const
114 {
115 #ifdef TEUCHOS_DEBUG
116 TEUCHOS_TEST_FOR_EXCEPTION(i > d,
117 std::out_of_range,
118 "Polynomial<CoeffT>::getCoefficient(i): " <<
119 "Error, coefficient i = " << i <<
120 " is not in range, degree = " << d << "." );
121 #endif
122 return coeff[i];
123 }
124
125 template <typename CoeffT>
126 void
setCoefficient(unsigned int i,const CoeffT & v)127 Teuchos::Polynomial<CoeffT>::setCoefficient(unsigned int i, const CoeffT& v)
128 {
129 #ifdef TEUCHOS_DEBUG
130 TEUCHOS_TEST_FOR_EXCEPTION(i > d,
131 std::out_of_range,
132 "Polynomial<CoeffT>::setCoefficient(i,v): " <<
133 "Error, coefficient i = " << i <<
134 " is not in range, degree = " << d << "." );
135 TEUCHOS_TEST_FOR_EXCEPTION(coeff[i] == Teuchos::null,
136 std::runtime_error,
137 "Polynomial<CoeffT>::setCoefficient(i,v): " <<
138 "Error, coefficient i = " << i << " is null!");
139 #endif
140 PolynomialTraits<CoeffT>::copy(v, coeff[i].get());
141 }
142
143 template <typename CoeffT>
144 void
setCoefficientPtr(unsigned int i,const Teuchos::RCP<CoeffT> & v)145 Teuchos::Polynomial<CoeffT>::setCoefficientPtr(
146 unsigned int i,
147 const Teuchos::RCP<CoeffT>& v)
148 {
149 #ifdef TEUCHOS_DEBUG
150 TEUCHOS_TEST_FOR_EXCEPTION(i > d,
151 std::out_of_range,
152 "Polynomial<CoeffT>::setCoefficientPtr(i,v): " <<
153 "Error, coefficient i = " << i <<
154 " is not in range, degree = " << d << "." );
155 #endif
156 coeff[i] = v;
157 }
158
159 template <typename CoeffT>
160 void
evaluate(typename Teuchos::Polynomial<CoeffT>::scalar_type & t,CoeffT * x,CoeffT * xdot) const161 Teuchos::Polynomial<CoeffT>::evaluate(
162 typename Teuchos::Polynomial<CoeffT>::scalar_type& t,
163 CoeffT* x, CoeffT* xdot) const
164 {
165 bool evaluate_xdot = (xdot != NULL);
166
167 #ifdef TEUCHOS_DEBUG
168 for (unsigned int i=0; i<=d; i++)
169 TEUCHOS_TEST_FOR_EXCEPTION(coeff[i] == Teuchos::null,
170 std::runtime_error,
171 "Polynomial<CoeffT>::evaluate(): " <<
172 "Error, coefficient i = " << i << " is null!");
173 #endif
174
175 // Initialize x, xdot with coeff[d]
176 PolynomialTraits<CoeffT>::copy(*coeff[d], x);
177 if (evaluate_xdot) {
178 if (d > 0)
179 PolynomialTraits<CoeffT>::copy(*coeff[d], xdot);
180 else
181 PolynomialTraits<CoeffT>::assign(
182 xdot,
183 Teuchos::ScalarTraits<scalar_type>::zero());
184 }
185
186 // If this is a degree 0 polynomial, we're done
187 if (d == 0)
188 return;
189
190 for (int k=d-1; k>=0; --k) {
191 // compute x = coeff[k] + t*x
192 PolynomialTraits<CoeffT>::update(x, *coeff[k], t);
193
194 // compute xdot = x + t*xdot
195 if (evaluate_xdot && k > 0)
196 PolynomialTraits<CoeffT>::update(xdot, *x, t);
197 }
198 }
199
200 #endif // TEUCHOS_VECTOR_POLYNOMIAL_HPP
201