← Back to branch summary

~ubuntu-branches/ubuntu/vivid/deal.ii/vivid

« back to all changes in this revision

Viewing changes to base/source/polynomial.cc

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV, Adam C. Powell, IV, Denis Barbier
  • Date: 2010-07-29 13:47:01 UTC
  • mfrom: (3.1.3 sid)
  • Revision ID: james.westby@ubuntu.com-20100729134701-akb8jb3stwge8tcm
Tags: 6.3.1-1
http://bugs.debian.org/562332
[ Adam C. Powell, IV ]
* Changed to source format 3.0 (quilt).
* Changed maintainer to debian-science with Adam Powell as uploader.
* Added source lintian overrides about Adam Powell's name.
* Added Vcs info on git repository.
* Bumped Standards-Version.
* Changed stamp-patch to patch target and fixed its application criterion.
* Moved make_dependencies and expand_instantiations to a versioned directory
  to avoid shlib package conflicts.

[ Denis Barbier ]
* New upstream release (closes: #562332).
  + Added libtbb support.
  + Forward-ported all patches.
* Updates for new PETSc version, including workaround for different versions
  of petsc and slepc.
* Add debian/watch.
* Update to debhelper 7.
* Added pdebuild patch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
base/source/polynomial.cc
1
1
//---------------------------------------------------------------------------
2
 
//      $Id: polynomial.cc 14038 2006-10-23 02:46:34Z bangerth $   
 
2
//      $Id: polynomial.cc 20941 2010-04-06 04:00:59Z bangerth $
3
3
//    Version: $Name$
4
4
//
5
 
//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
 
5
//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010 by the deal.II authors
6
6
//
7
7
//    This file is subject to QPL and may not be  distributed
8
8
//    without copyright and license information. Please refer
31
31
// to be a problem since we only need it on very rare occasions. if
32
32
// someone finds this is a bottleneck, feel free to replace it by a
33
33
// more fine-grained solution
34
 
namespace 
 
34
namespace
35
35
{
36
36
  Threads::ThreadMutex coefficients_lock;
37
37
}
53
53
 
54
54
 
55
55
  template <typename number>
56
 
  number
57
 
  Polynomial<number>::value (const number x) const
58
 
  {
59
 
    Assert (coefficients.size() > 0, ExcEmptyObject());
60
 
    const unsigned int m=coefficients.size();
61
 
 
62
 
                                     // Horner scheme
63
 
    number value = coefficients.back();
64
 
    for (int k=m-2; k>=0; --k)
65
 
      value = value*x + coefficients[k];
66
 
 
67
 
    return value;
68
 
  }
 
56
  Polynomial<number>::Polynomial (const unsigned int n)
 
57
                  :
 
58
                  coefficients(n+1, 0.)
 
59
  {}
69
60
 
70
61
 
71
62
 
77
68
    Assert (coefficients.size() > 0, ExcEmptyObject());
78
69
    Assert (values.size() > 0, ExcZero());
79
70
    const unsigned int values_size=values.size();
80
 
  
81
 
  
 
71
 
 
72
 
82
73
                                     // if we only need the value, then
83
74
                                     // call the other function since
84
75
                                     // that is significantly faster
133
124
      {
134
125
        *c *= f;
135
126
        f *= factor;
136
 
      }  
 
127
      }
137
128
  }
138
129
 
139
130
 
169
160
    return *this;
170
161
  }
171
162
 
172
 
  
 
163
 
173
164
  template <typename number>
174
165
  Polynomial<number>&
175
166
  Polynomial<number>::operator *= (const Polynomial<number>& p)
178
169
    unsigned int new_degree = this->degree() + p.degree();
179
170
 
180
171
    std::vector<number> new_coefficients(new_degree+1, 0.);
181
 
    
 
172
 
182
173
    for (unsigned int i=0; i<p.coefficients.size(); ++i)
183
174
      for (unsigned int j=0; j<this->coefficients.size(); ++j)
184
175
      new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i];
185
176
    this->coefficients = new_coefficients;
186
 
    
 
177
 
187
178
    return *this;
188
179
  }
189
180
 
190
 
  
 
181
 
191
182
  template <typename number>
192
183
  Polynomial<number>&
193
184
  Polynomial<number>::operator += (const Polynomial<number>& p)
203
194
    return *this;
204
195
  }
205
196
 
206
 
  
 
197
 
207
198
  template <typename number>
208
199
  Polynomial<number>&
209
200
  Polynomial<number>::operator -= (const Polynomial<number>& p)
219
210
    return *this;
220
211
  }
221
212
 
222
 
  
 
213
 
223
214
  template <typename number>
224
215
  template <typename number2>
225
216
  void
226
217
  Polynomial<number>::shift(std::vector<number>& coefficients,
227
218
                            const number2 offset)
228
 
  {  
 
219
  {
229
220
#ifdef DEAL_II_LONG_DOUBLE_LOOP_BUG
230
221
    AssertThrow (false,
231
222
                 ExcMessage("Sorry, but the compiler you are using has a bug that disallows "
237
228
                                     // args. note that this code is
238
229
                                     // actually unreachable
239
230
    coefficients[0] = offset;
240
 
#else  
 
231
#else
241
232
                                     // Copy coefficients to a vector of
242
233
                                     // accuracy given by the argument
243
234
    std::vector<number2> new_coefficients(coefficients.begin(),
244
235
                                          coefficients.end());
245
 
  
 
236
 
246
237
                                     // Traverse all coefficients from
247
238
                                     // c_1. c_0 will be modified by
248
239
                                     // higher degrees, only.
259
250
                                         // needed and computed
260
251
                                         // successively.
261
252
        number2 offset_power = offset;
262
 
      
 
253
 
263
254
                                         // Compute (x+offset)^d
264
255
                                         // and modify all values c_k
265
256
                                         // with k<d.
300
291
  }
301
292
 
302
293
 
303
 
  
 
294
 
304
295
  template <typename number>
305
296
  Polynomial<number>
306
297
  Polynomial<number>::derivative () const
314
305
 
315
306
    return Polynomial<number> (newcoefficients);
316
307
  }
317
 
  
 
308
 
318
309
 
319
310
  template <typename number>
320
311
  Polynomial<number>
327
318
 
328
319
    return Polynomial<number> (newcoefficients);
329
320
  }
330
 
  
 
321
 
331
322
 
332
323
  template <typename number>
333
324
  void
352
343
    result[n] = coefficient;
353
344
    return result;
354
345
  }
355
 
  
356
 
  
 
346
 
 
347
 
357
348
  template <typename number>
358
349
  Monomial<number>::Monomial (unsigned int n,
359
350
                              double coefficient)
360
351
                  : Polynomial<number>(make_vector(n, coefficient))
361
352
  {}
362
 
  
363
 
  
 
353
 
 
354
 
364
355
  template <typename number>
365
356
  std::vector<Polynomial<number> >
366
357
  Monomial<number>::generate_complete_basis (const unsigned int degree)
375
366
 
376
367
  LagrangeEquidistant::LagrangeEquidistant (const unsigned int n,
377
368
                                            const unsigned int support_point)
378
 
                  :
379
 
                  Polynomial<double>(compute_coefficients(n,support_point))
380
 
  {}
381
 
 
382
 
 
383
 
 
384
 
  std::vector<double> 
 
369
  {
 
370
    if (n <= 10)
 
371
      {
 
372
        this->coefficients.resize(n+1);
 
373
        compute_coefficients(n, support_point, this->coefficients);
 
374
      }
 
375
    else
 
376
      {
 
377
                                         // We have precomputed tables
 
378
                                         // up to degree 10. For
 
379
                                         // higher order, we have to
 
380
                                         // compute by hand.
 
381
 
 
382
                                         // Start with the constant one
 
383
        this->coefficients.resize(1);
 
384
        this->coefficients[0] = 1.;
 
385
 
 
386
                                         // Then compute the Lagrange
 
387
                                         // polynomial as the product
 
388
                                         // of linear factors
 
389
        std::vector<double> two (2, 1.);
 
390
 
 
391
        for (unsigned int k=0;k<=n;++k)
 
392
          {
 
393
            if (k != support_point)
 
394
              {
 
395
                two[0] = -1.*k/n;
 
396
                Polynomial<double> factor(two);
 
397
                factor.scale(1.*n/(support_point - k));
 
398
                (*this) *= factor;
 
399
              }
 
400
          }
 
401
      }
 
402
  }
 
403
 
 
404
 
 
405
  void
385
406
  LagrangeEquidistant::compute_coefficients (const unsigned int n,
386
 
                                             const unsigned int support_point)
 
407
                                             const unsigned int support_point,
 
408
                                             std::vector<double>& a)
387
409
  {
388
 
    std::vector<double> a (n+1);
389
410
    Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));
390
411
 
391
412
    unsigned int n_functions=n+1;
392
413
    Assert(support_point<n_functions,
393
414
           ExcIndexRange(support_point, 0, n_functions));
394
415
    double const *x=0;
395
 
  
 
416
 
396
417
    switch (n)
397
418
      {
398
419
        case 1:
403
424
                  0.0, 1.0
404
425
            };
405
426
          x=&x1[0];
406
 
          break;        
 
427
          break;
407
428
        }
408
429
        case 2:
409
430
        {
439
460
                  0.0, -1.0, 22.0/3.0, -16.0, 32.0/3.0
440
461
            };
441
462
          x=&x4[0];
442
 
          break;        
 
463
          break;
443
464
        }
444
465
        case 5:
445
466
        {
592
613
          break;
593
614
        }
594
615
        default:
595
 
              Assert(false, ExcNotImplemented());
 
616
              Assert(false, ExcInternalError())
596
617
      }
597
618
 
598
619
    Assert(x!=0, ExcInternalError());
599
620
    for (unsigned int i=0; i<n_functions; ++i)
600
621
      a[i]=x[support_point*n_functions+i];
601
 
  
602
 
    return a;
603
622
  }
604
623
 
605
624
 
625
644
 
626
645
//----------------------------------------------------------------------//
627
646
 
628
 
  
 
647
 
629
648
  std::vector<Polynomial<double> >
630
649
  Lagrange::generate_complete_basis (const std::vector<Point<1> >& points)
631
650
  {
638
657
    std::vector<double> linear(2, 1.);
639
658
                                     // We start with a constant polynomial
640
659
    std::vector<double> one(1, 1.);
641
 
    
 
660
 
642
661
    for (unsigned int i=0;i<p.size();++i)
643
662
      {
644
663
                                         // Construct interpolation formula
675
694
              }
676
695
          }
677
696
      }
678
 
    
 
697
 
679
698
    return p;
680
699
  }
681
 
  
 
700
 
682
701
 
683
702
// ------------------ class Legendre --------------- //
684
703
 
685
704
 
686
 
//TODO:[?] This class leaks memory, but only at the very end of a program.
687
 
// Since it expands the Legendre<number>::coefficients array, the elements
688
 
// of this static variable are not destroyed at the end of the program
689
 
// run. While this is not a problem (since the returned memory could
690
 
// not be used anyway then), it is a little confusing when looking at
691
 
// a memory checker such as "purify". Maybe, this should be handled somehow
692
 
// to avoid this confusion in future.
693
 
 
694
705
// Reserve space for polynomials up to degree 19. Should be sufficient
695
706
// for the start.
696
 
  std::vector<const std::vector<double> *>
697
 
  Legendre::recursive_coefficients(20,
698
 
                                   static_cast<const std::vector<double>*>(0));
699
 
  std::vector<const std::vector<double> *>
700
 
  Legendre::shifted_coefficients(20,
701
 
                                 static_cast<const std::vector<double>*>(0));
 
707
  std::vector<std_cxx1x::shared_ptr<const std::vector<double> > >
 
708
  Legendre::recursive_coefficients(20);
 
709
  std::vector<std_cxx1x::shared_ptr<const std::vector<double> > >
 
710
  Legendre::shifted_coefficients(20);
702
711
 
703
712
 
704
713
  Legendre::Legendre (const unsigned int k)
707
716
  {}
708
717
 
709
718
 
710
 
  
 
719
 
711
720
  void
712
721
  Legendre::compute_coefficients (const unsigned int k_)
713
722
  {
721
730
#else
722
731
    typedef long double SHIFT_TYPE;
723
732
#endif
724
 
    
 
733
 
725
734
    unsigned int k = k_;
726
735
 
727
736
                                     // first make sure that no other
738
747
                                     // already exist?
739
748
    if ((recursive_coefficients.size() < k+1) ||
740
749
        ((recursive_coefficients.size() >= k+1) &&
741
 
         (recursive_coefficients[k] == 0)))
 
750
         (recursive_coefficients[k] ==
 
751
          std_cxx1x::shared_ptr<const std::vector<double> >())))
742
752
                                       // no, then generate the
743
753
                                       // respective coefficients
744
754
      {
745
 
        recursive_coefficients.resize (k+1, 0);
746
 
      
 
755
        recursive_coefficients.resize (k+1);
 
756
 
747
757
        if (k<=1)
748
758
          {
749
759
                                             // create coefficients
762
772
            (*c1)[1] = 1.;
763
773
 
764
774
                                             // now make these arrays
765
 
                                             // const
766
 
            recursive_coefficients[0] = c0;
767
 
            recursive_coefficients[1] = c1;
 
775
                                             // const. use shared_ptr for
 
776
                                             // recursive_coefficients because
 
777
                                             // that avoids a memory leak that
 
778
                                             // would appear if we used plain
 
779
                                             // pointers.
 
780
            recursive_coefficients[0] =
 
781
              std_cxx1x::shared_ptr<const std::vector<double> >(c0);
 
782
            recursive_coefficients[1] =
 
783
              std_cxx1x::shared_ptr<const std::vector<double> >(c1);
 
784
 
768
785
                                             // Compute polynomials
769
786
                                             // orthogonal on [0,1]
770
787
            c0 = new std::vector<double>(*c0);
771
788
            c1 = new std::vector<double>(*c1);
772
 
          
 
789
 
773
790
            Polynomial<double>::shift<SHIFT_TYPE> (*c0, -1.);
774
791
            Polynomial<double>::scale(*c0, 2.);
775
792
            Polynomial<double>::shift<SHIFT_TYPE> (*c1, -1.);
776
793
            Polynomial<double>::scale(*c1, 2.);
777
794
            Polynomial<double>::multiply(*c1, std::sqrt(3.));
778
 
            shifted_coefficients[0]=c0;
779
 
            shifted_coefficients[1]=c1;
 
795
            shifted_coefficients[0]=std_cxx1x::shared_ptr<const std::vector<double> >(c0);
 
796
            shifted_coefficients[1]=std_cxx1x::shared_ptr<const std::vector<double> >(c1);
780
797
          }
781
798
        else
782
799
          {
793
810
            coefficients_lock.acquire ();
794
811
 
795
812
            std::vector<double> *ck = new std::vector<double>(k+1);
796
 
          
 
813
 
797
814
            const double a = 1./(k);
798
815
            const double b = a*(2*k-1);
799
816
            const double c = a*(k-1);
800
 
          
 
817
 
801
818
            (*ck)[k]   = b*(*recursive_coefficients[k-1])[k-1];
802
819
            (*ck)[k-1] = b*(*recursive_coefficients[k-1])[k-2];
803
820
            for (unsigned int i=1 ; i<= k-2 ; ++i)
810
827
                                             // created vector to the
811
828
                                             // const pointer in the
812
829
                                             // coefficients array
813
 
            recursive_coefficients[k] = ck;
 
830
            recursive_coefficients[k] =
 
831
              std_cxx1x::shared_ptr<const std::vector<double> >(ck);
814
832
                                             // and compute the
815
833
                                             // coefficients for [0,1]
816
834
            ck = new std::vector<double>(*ck);
817
835
            Polynomial<double>::shift<SHIFT_TYPE> (*ck, -1.);
818
836
            Polynomial<double>::scale(*ck, 2.);
819
837
            Polynomial<double>::multiply(*ck, std::sqrt(2.*k+1.));
820
 
            shifted_coefficients[k] = ck;
 
838
            shifted_coefficients[k] =
 
839
              std_cxx1x::shared_ptr<const std::vector<double> >(ck);
821
840
          };
822
841
      };
823
842
  }
883
902
                                     // until we quit this function
884
903
    Threads::ThreadMutex::ScopedLock lock(coefficients_lock);
885
904
 
886
 
                                     // The first 2 coefficients 
 
905
                                     // The first 2 coefficients
887
906
                                     // are hard-coded
888
907
    if (k==0)
889
908
      k=1;
890
909
                                     // check: does the information
891
910
                                     // already exist?
892
911
    if (  (recursive_coefficients.size() < k+1) ||
893
 
          ((recursive_coefficients.size() >= k+1) && 
 
912
          ((recursive_coefficients.size() >= k+1) &&
894
913
           (recursive_coefficients[k] == 0)) )
895
914
                                           // no, then generate the
896
915
                                           // respective coefficients
897
916
      {
898
917
        recursive_coefficients.resize (k+1, 0);
899
 
      
 
918
 
900
919
        if (k<=1)
901
920
          {
902
921
                                             // create coefficients
933
952
            (*c2)[0] =   0.*a;
934
953
            (*c2)[1] =  -4.*a;
935
954
            (*c2)[2] =   4.*a;
936
 
            
 
955
 
937
956
            recursive_coefficients[2] = c2;
938
957
          }
939
958
        else
951
970
            coefficients_lock.acquire ();
952
971
 
953
972
            std::vector<double> *ck = new std::vector<double>(k+1);
954
 
           
 
973
 
955
974
            const double a = 1.; //1./(2.*k);
956
975
 
957
976
            (*ck)[0] = - a*(*recursive_coefficients[k-1])[0];
958
 
          
 
977
 
959
978
            for (unsigned int i=1; i<=k-1; ++i)
960
979
                (*ck)[i] = a*( 2.*(*recursive_coefficients[k-1])[i-1]
961
980
                               - (*recursive_coefficients[k-1])[i] );
962
 
          
 
981
 
963
982
            (*ck)[k] = a*2.*(*recursive_coefficients[k-1])[k-1];
964
983
                                          // for even degrees, we need
965
984
                                          // to add a multiple of
972
991
 
973
992
                (*ck)[1] += b*(*recursive_coefficients[2])[1];
974
993
                (*ck)[2] += b*(*recursive_coefficients[2])[2];
975
 
              }   
 
994
              }
976
995
                                             // finally assign the newly
977
996
                                             // created vector to the
978
997
                                             // const pointer in the