]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make compile with gcc 3.3.x. It has trouble finding the proper
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 02:44:35 +0000 (02:44 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 02:44:35 +0000 (02:44 +0000)
overload of std::sqrt in some cases.

git-svn-id: https://svn.dealii.org/trunk@21971 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/polynomial.cc

index e82e048956444dc67c0c1e737e739c2772d426ad..a84e6379b3767afcb650711d37694d15cc0a78af 100644 (file)
@@ -1,8 +1,8 @@
 //---------------------------------------------------------------------------
-//      $Id$   
+//      $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -31,7 +31,7 @@ DEAL_II_NAMESPACE_OPEN
 // to be a problem since we only need it on very rare occasions. if
 // someone finds this is a bottleneck, feel free to replace it by a
 // more fine-grained solution
-namespace 
+namespace
 {
   Threads::ThreadMutex coefficients_lock;
 }
@@ -60,8 +60,8 @@ namespace Polynomials
     Assert (coefficients.size() > 0, ExcEmptyObject());
     Assert (values.size() > 0, ExcZero());
     const unsigned int values_size=values.size();
-  
-  
+
+
                                      // if we only need the value, then
                                      // call the other function since
                                      // that is significantly faster
@@ -116,7 +116,7 @@ namespace Polynomials
       {
         *c *= f;
         f *= factor;
-      }  
+      }
   }
 
 
@@ -152,7 +152,7 @@ namespace Polynomials
     return *this;
   }
 
-  
+
   template <typename number>
   Polynomial<number>&
   Polynomial<number>::operator *= (const Polynomial<number>& p)
@@ -161,16 +161,16 @@ namespace Polynomials
     unsigned int new_degree = this->degree() + p.degree();
 
     std::vector<number> new_coefficients(new_degree+1, 0.);
-    
+
     for (unsigned int i=0; i<p.coefficients.size(); ++i)
       for (unsigned int j=0; j<this->coefficients.size(); ++j)
       new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i];
     this->coefficients = new_coefficients;
-    
+
     return *this;
   }
 
-  
+
   template <typename number>
   Polynomial<number>&
   Polynomial<number>::operator += (const Polynomial<number>& p)
@@ -186,7 +186,7 @@ namespace Polynomials
     return *this;
   }
 
-  
+
   template <typename number>
   Polynomial<number>&
   Polynomial<number>::operator -= (const Polynomial<number>& p)
@@ -202,13 +202,13 @@ namespace Polynomials
     return *this;
   }
 
-  
+
   template <typename number>
   template <typename number2>
   void
   Polynomial<number>::shift(std::vector<number>& coefficients,
                             const number2 offset)
-  {  
+  {
 #ifdef DEAL_II_LONG_DOUBLE_LOOP_BUG
     AssertThrow (false,
                  ExcMessage("Sorry, but the compiler you are using has a bug that disallows "
@@ -220,12 +220,12 @@ namespace Polynomials
                                      // args. note that this code is
                                      // actually unreachable
     coefficients[0] = offset;
-#else  
+#else
                                      // Copy coefficients to a vector of
                                      // accuracy given by the argument
     std::vector<number2> new_coefficients(coefficients.begin(),
                                           coefficients.end());
-  
+
                                      // Traverse all coefficients from
                                      // c_1. c_0 will be modified by
                                      // higher degrees, only.
@@ -242,7 +242,7 @@ namespace Polynomials
                                          // needed and computed
                                          // successively.
         number2 offset_power = offset;
-      
+
                                          // Compute (x+offset)^d
                                          // and modify all values c_k
                                          // with k<d.
@@ -283,7 +283,7 @@ namespace Polynomials
   }
 
 
-  
+
   template <typename number>
   Polynomial<number>
   Polynomial<number>::derivative () const
@@ -297,7 +297,7 @@ namespace Polynomials
 
     return Polynomial<number> (newcoefficients);
   }
-  
+
 
   template <typename number>
   Polynomial<number>
@@ -310,7 +310,7 @@ namespace Polynomials
 
     return Polynomial<number> (newcoefficients);
   }
-  
+
 
   template <typename number>
   void
@@ -335,15 +335,15 @@ namespace Polynomials
     result[n] = coefficient;
     return result;
   }
-  
-  
+
+
   template <typename number>
   Monomial<number>::Monomial (unsigned int n,
                              double coefficient)
                  : Polynomial<number>(make_vector(n, coefficient))
   {}
-  
-  
+
+
   template <typename number>
   std::vector<Polynomial<number> >
   Monomial<number>::generate_complete_basis (const unsigned int degree)
@@ -364,7 +364,7 @@ namespace Polynomials
 
 
 
-  std::vector<double> 
+  std::vector<double>
   LagrangeEquidistant::compute_coefficients (const unsigned int n,
                                              const unsigned int support_point)
   {
@@ -375,7 +375,7 @@ namespace Polynomials
     Assert(support_point<n_functions,
            ExcIndexRange(support_point, 0, n_functions));
     double const *x=0;
-  
+
     switch (n)
       {
         case 1:
@@ -386,7 +386,7 @@ namespace Polynomials
                   0.0, 1.0
             };
           x=&x1[0];
-          break;       
+          break;
         }
         case 2:
         {
@@ -422,7 +422,7 @@ namespace Polynomials
                   0.0, -1.0, 22.0/3.0, -16.0, 32.0/3.0
             };
           x=&x4[0];
-          break;       
+          break;
         }
         case 5:
         {
@@ -581,7 +581,7 @@ namespace Polynomials
     Assert(x!=0, ExcInternalError());
     for (unsigned int i=0; i<n_functions; ++i)
       a[i]=x[support_point*n_functions+i];
-  
+
     return a;
   }
 
@@ -608,7 +608,7 @@ namespace Polynomials
 
 //----------------------------------------------------------------------//
 
-  
+
   std::vector<Polynomial<double> >
   Lagrange::generate_complete_basis (const std::vector<Point<1> >& points)
   {
@@ -621,7 +621,7 @@ namespace Polynomials
     std::vector<double> linear(2, 1.);
                                     // We start with a constant polynomial
     std::vector<double> one(1, 1.);
-    
+
     for (unsigned int i=0;i<p.size();++i)
       {
                                         // Construct interpolation formula
@@ -658,10 +658,10 @@ namespace Polynomials
              }
          }
       }
-    
+
     return p;
   }
-  
+
 
 // ------------------ class Legendre --------------- //
 
@@ -690,7 +690,7 @@ namespace Polynomials
   {}
 
 
-  
+
   void
   Legendre::compute_coefficients (const unsigned int k_)
   {
@@ -704,7 +704,7 @@ namespace Polynomials
 #else
     typedef long double SHIFT_TYPE;
 #endif
-    
+
     unsigned int k = k_;
 
                                      // first make sure that no other
@@ -726,7 +726,7 @@ namespace Polynomials
                                        // respective coefficients
       {
         recursive_coefficients.resize (k+1, 0);
-      
+
         if (k<=1)
           {
                                              // create coefficients
@@ -752,7 +752,7 @@ namespace Polynomials
                                              // orthogonal on [0,1]
             c0 = new std::vector<double>(*c0);
             c1 = new std::vector<double>(*c1);
-         
+
             Polynomial<double>::shift<SHIFT_TYPE> (*c0, -1.);
             Polynomial<double>::scale(*c0, 2.);
             Polynomial<double>::shift<SHIFT_TYPE> (*c1, -1.);
@@ -776,11 +776,11 @@ namespace Polynomials
             coefficients_lock.acquire ();
 
             std::vector<double> *ck = new std::vector<double>(k+1);
-         
+
             const double a = 1./(k);
             const double b = a*(2*k-1);
             const double c = a*(k-1);
-         
+
             (*ck)[k]   = b*(*recursive_coefficients[k-1])[k-1];
             (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2];
             for (unsigned int i=1 ; i<= k-2 ; ++i)
@@ -862,8 +862,8 @@ std::vector<double> Lobatto::compute_coefficients (const unsigned int p) {
          std::vector<double> coefficients (3);
 
          coefficients[0] = 0.0;
-         coefficients[1] = -1.0 * std::sqrt (3);
-         coefficients[2] = std::sqrt (3);
+         coefficients[1] = -1.0 * std::sqrt (3.);
+         coefficients[2] = std::sqrt (3.);
          return coefficients;
       }
 
@@ -872,8 +872,8 @@ std::vector<double> Lobatto::compute_coefficients (const unsigned int p) {
          std::vector<double> legendre_coefficients_tmp1 (p);
          std::vector<double> legendre_coefficients_tmp2 (p - 1);
 
-         coefficients[0] = -1.0 * std::sqrt (3);
-         coefficients[1] = 2.0 * std::sqrt (3);
+         coefficients[0] = -1.0 * std::sqrt (3.);
+         coefficients[1] = 2.0 * std::sqrt (3.);
          legendre_coefficients_tmp1[0] = 1.0;
 
          for (unsigned int i = 2; i < p; ++i) {
@@ -883,13 +883,13 @@ std::vector<double> Lobatto::compute_coefficients (const unsigned int p) {
             for (unsigned int j = 0; j < i; ++j)
                legendre_coefficients_tmp1[j] = coefficients[j];
 
-            coefficients[0] = std::sqrt (2 * i + 1) * ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] / std::sqrt (2 * i - 1) + (1.0 - i) * legendre_coefficients_tmp2[0] / std::sqrt (2 * i - 3)) / i;
+            coefficients[0] = std::sqrt (2 * i + 1.) * ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] / std::sqrt (2 * i - 1.) + (1.0 - i) * legendre_coefficients_tmp2[0] / std::sqrt (2 * i - 3.)) / i;
 
             for (unsigned int j = 1; j < i - 1; ++j)
-               coefficients[j] = std::sqrt (2 * i + 1) * (std::sqrt (2 * i - 1) * (2.0 * legendre_coefficients_tmp1[j - 1] - legendre_coefficients_tmp1[j]) + (1.0 - i) * legendre_coefficients_tmp2[j] / std::sqrt (2 * i - 3)) / i;
+               coefficients[j] = std::sqrt (2 * i + 1.) * (std::sqrt (2 * i - 1.) * (2.0 * legendre_coefficients_tmp1[j - 1] - legendre_coefficients_tmp1[j]) + (1.0 - i) * legendre_coefficients_tmp2[j] / std::sqrt (2 * i - 3.)) / i;
 
-            coefficients[i - 1] = std::sqrt (4 * i * i - 1) * (2.0 * legendre_coefficients_tmp1[i - 2] - legendre_coefficients_tmp1[i - 1]) / i;
-            coefficients[i] = 2.0 * std::sqrt (4 * i * i - 1) * legendre_coefficients_tmp1[i - 1] / i;
+            coefficients[i - 1] = std::sqrt (4 * i * i - 1.) * (2.0 * legendre_coefficients_tmp1[i - 2] - legendre_coefficients_tmp1[i - 1]) / i;
+            coefficients[i] = 2.0 * std::sqrt (4 * i * i - 1.) * legendre_coefficients_tmp1[i - 1] / i;
          }
 
          for (int i = p; i > 0; --i)
@@ -941,20 +941,20 @@ std::vector<Polynomial<double> > Lobatto::generate_complete_basis (const unsigne
                                      // until we quit this function
     Threads::ThreadMutex::ScopedLock lock(coefficients_lock);
 
-                                     // The first 2 coefficients 
+                                     // The first 2 coefficients
                                      // are hard-coded
     if (k==0)
       k=1;
                                      // check: does the information
                                      // already exist?
     if (  (recursive_coefficients.size() < k+1) ||
-         ((recursive_coefficients.size() >= k+1) && 
+         ((recursive_coefficients.size() >= k+1) &&
            (recursive_coefficients[k] == 0)) )
                                           // no, then generate the
                                           // respective coefficients
       {
        recursive_coefficients.resize (k+1, 0);
-      
+
        if (k<=1)
          {
                                              // create coefficients
@@ -991,7 +991,7 @@ std::vector<Polynomial<double> > Lobatto::generate_complete_basis (const unsigne
            (*c2)[0] =   0.*a;
            (*c2)[1] =  -4.*a;
            (*c2)[2] =   4.*a;
-           
+
            recursive_coefficients[2] = c2;
          }
        else
@@ -1009,15 +1009,15 @@ std::vector<Polynomial<double> > Lobatto::generate_complete_basis (const unsigne
            coefficients_lock.acquire ();
 
            std::vector<double> *ck = new std::vector<double>(k+1);
-          
+
            const double a = 1.; //1./(2.*k);
 
            (*ck)[0] = - a*(*recursive_coefficients[k-1])[0];
-         
+
            for (unsigned int i=1; i<=k-1; ++i)
                (*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1]
                               - (*recursive_coefficients[k-1])[i] );
-         
+
            (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1];
                                          // for even degrees, we need
                                          // to add a multiple of
@@ -1030,7 +1030,7 @@ std::vector<Polynomial<double> > Lobatto::generate_complete_basis (const unsigne
 
                (*ck)[1] += b*(*recursive_coefficients[2])[1];
                (*ck)[2] += b*(*recursive_coefficients[2])[2];
-             }   
+             }
                                              // finally assign the newly
                                              // created vector to the
                                              // const pointer in the

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.