Teuchos_Polynomial.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef TEUCHOS_POLYNOMIAL_HPP
00030 #define TEUCHOS_POLYNOMIAL_HPP
00031 
00032 #include "Teuchos_PolynomialDecl.hpp"
00033 #include "Teuchos_ScalarTraits.hpp"
00034 
00035 template <typename CoeffT>
00036 Teuchos::Polynomial<CoeffT>::Polynomial(unsigned int deg,
00037           const CoeffT& cloneCoeff, 
00038           unsigned int reserve) :
00039   d(deg)
00040 {
00041   if (reserve > d)
00042     sz = reserve+1;
00043   else
00044     sz = d+1;
00045 
00046   coeff.resize(sz);
00047   for (unsigned int i=0; i<sz; i++)
00048     coeff[i] = PolynomialTraits<CoeffT>::clone(cloneCoeff);
00049 }
00050 
00051 template <typename CoeffT>
00052 Teuchos::Polynomial<CoeffT>::Polynomial(unsigned int deg, 
00053           unsigned int reserve) :
00054   d(deg)
00055 {
00056   if (reserve > d)
00057     sz = reserve+1;
00058   else
00059     sz = d+1;
00060 
00061   coeff.resize(sz);
00062 }
00063 
00064 template <typename CoeffT>
00065 Teuchos::Polynomial<CoeffT>::~Polynomial()
00066 {
00067 }
00068 
00069 template <typename CoeffT>
00070 void 
00071 Teuchos::Polynomial<CoeffT>::setDegree(unsigned int deg) 
00072 {
00073   d = deg;
00074   if (d+1 > sz) {
00075     coeff.resize(d+1);
00076     if (coeff[0] != Teuchos::null) {
00077       for (unsigned int i=sz; i<d+1; i++)
00078   coeff[i] = PolynomialTraits<CoeffT>::clone(*coeff[0]);
00079     }
00080     sz = d+1;
00081   }
00082 }
00083 
00084 template <typename CoeffT>
00085 Teuchos::RCP<CoeffT>
00086 Teuchos::Polynomial<CoeffT>::getCoefficient(unsigned int i) 
00087 {
00088 #ifdef TEUCHOS_DEBUG
00089   TEST_FOR_EXCEPTION(i > d, 
00090          std::out_of_range,
00091          "Polynomial<CoeffT>::getCoefficient(i): " << 
00092          "Error, coefficient i = " << i << 
00093          " is not in range, degree = " << d << "." );
00094 #endif
00095   return coeff[i];
00096 }
00097 
00098 template <typename CoeffT>
00099 Teuchos::RCP<const CoeffT>
00100 Teuchos::Polynomial<CoeffT>::getCoefficient(unsigned int i) const
00101 {
00102 #ifdef TEUCHOS_DEBUG
00103   TEST_FOR_EXCEPTION(i > d, 
00104          std::out_of_range,
00105          "Polynomial<CoeffT>::getCoefficient(i): " << 
00106          "Error, coefficient i = " << i << 
00107          " is not in range, degree = " << d << "." );
00108 #endif
00109   return coeff[i];
00110 }
00111 
00112 template <typename CoeffT>
00113 void
00114 Teuchos::Polynomial<CoeffT>::setCoefficient(unsigned int i, const CoeffT& v) 
00115 {
00116 #ifdef TEUCHOS_DEBUG
00117   TEST_FOR_EXCEPTION(i > d, 
00118          std::out_of_range,
00119          "Polynomial<CoeffT>::setCoefficient(i,v): " << 
00120          "Error, coefficient i = " << i << 
00121          " is not in range, degree = " << d << "." );
00122   TEST_FOR_EXCEPTION(coeff[i] == Teuchos::null, 
00123          std::runtime_error,
00124          "Polynomial<CoeffT>::setCoefficient(i,v): " << 
00125          "Error, coefficient i = " << i << " is null!");
00126 #endif
00127   PolynomialTraits<CoeffT>::copy(v, coeff[i].get());
00128 }
00129 
00130 template <typename CoeffT>
00131 void
00132 Teuchos::Polynomial<CoeffT>::setCoefficientPtr(
00133                             unsigned int i,
00134                             const Teuchos::RCP<CoeffT>& v)
00135 {
00136 #ifdef TEUCHOS_DEBUG
00137   TEST_FOR_EXCEPTION(i > d, 
00138          std::out_of_range,
00139          "Polynomial<CoeffT>::setCoefficientPtr(i,v): " << 
00140          "Error, coefficient i = " << i << 
00141          " is not in range, degree = " << d << "." );
00142 #endif
00143   coeff[i] = v;
00144 }
00145 
00146 template <typename CoeffT>
00147 void
00148 Teuchos::Polynomial<CoeffT>::evaluate(
00149         typename Teuchos::Polynomial<CoeffT>::scalar_type t, 
00150         CoeffT* x, CoeffT* xdot) const
00151 {
00152   bool evaluate_xdot = (xdot != NULL);
00153 
00154 #ifdef TEUCHOS_DEBUG
00155   for (unsigned int i=0; i<=d; i++)
00156     TEST_FOR_EXCEPTION(coeff[i] == Teuchos::null, 
00157            std::runtime_error,
00158            "Polynomial<CoeffT>::evaluate(): " << 
00159            "Error, coefficient i = " << i << " is null!");
00160 #endif
00161 
00162   // Initialize x, xdot with coeff[d]
00163   PolynomialTraits<CoeffT>::copy(*coeff[d], x);
00164   if (evaluate_xdot) {
00165     if (d > 0)
00166       PolynomialTraits<CoeffT>::copy(*coeff[d], xdot);
00167     else
00168       PolynomialTraits<CoeffT>::assign(
00169          xdot, 
00170          Teuchos::ScalarTraits<scalar_type>::zero());
00171   }
00172 
00173   // If this is a degree 0 polynomial, we're done
00174   if (d == 0)
00175     return;
00176 
00177   for (int k=d-1; k>=0; --k) {
00178     // compute x = coeff[k] + t*x
00179     PolynomialTraits<CoeffT>::update(x, *coeff[k], t);
00180 
00181     // compute xdot = x + t*xdot
00182     if (evaluate_xdot && k > 0)
00183       PolynomialTraits<CoeffT>::update(xdot, *x, t);
00184   }
00185 }
00186 
00187 #endif  // TEUCHOS_VECTOR_POLYNOMIAL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 09:57:30 2011 for Teuchos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3