Rheolef
7.1
an efficient C++ finite element environment
gauss_jacobi.icc
Go to the documentation of this file.
1
#include "rheolef/compiler.h"
22
#include "rheolef/gamma.h"
23
#include "rheolef/jacobi.h"
24
#include "rheolef/jacobi_roots.h"
25
#include <iterator>
26
namespace
rheolef
{
27
template
<
class
Size,
class
OutputIterator1,
class
OutputIterator2>
28
void
29
gauss_jacobi
(Size R,
30
typename
std::iterator_traits<OutputIterator1>::value_type
alpha
,
31
typename
std::iterator_traits<OutputIterator1>::value_type
beta
,
32
OutputIterator1 zeta, OutputIterator2 omega)
33
{
34
typedef
typename
std::iterator_traits<OutputIterator1>::value_type
T
;
35
T
num =
pow
(
T
(2.),
alpha
+
beta
+3)/sqr(
alpha
+
beta
+
T
(1.*R)+1);
36
if
(
alpha
== floor(
alpha
) &&
beta
== floor(
beta
))
37
for
(Size k = 1; k <= size_t(static_cast<int>(
beta
)); k++)
38
num *= (
T
(1.*R)+
T
(1.*k))/(
alpha
+
T
(1.*R)+
T
(1.*k));
39
else
40
num *= (
my_gamma
(
alpha
+
T
(1.*R)+1)/
my_gamma
(
alpha
+
beta
+
T
(1.*R)+1))
41
*(
my_gamma
(
beta
+
T
(1.*R)+1)/
my_gamma
(
T
(1.*R)+1));
42
jacobi_roots
(R,
alpha
,
beta
, zeta);
43
jacobi<T>
P (R-1,
alpha
+1,
beta
+1);
44
for
(Size r = 0; r < R; r++)
45
omega[r] = num/((1-sqr(zeta[r]))*sqr(P(zeta[r])));
46
}
47
}
// namespace rheolef
rheolef::my_gamma
T my_gamma(const T &x)
Definition:
gamma.icc:25
rheolef::gauss_jacobi
void gauss_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)
Definition:
gauss_jacobi.icc:29
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition:
bdf.icc:28
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition:
space_mult.h:120
rheolef
This file is part of Rheolef.
Definition:
compiler_eigen.h:37
rheolef::jacobi
Definition:
jacobi.icc:24
rheolef::jacobi_roots
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
Definition:
jacobi_roots.icc:26
rk::beta
Float beta[][pmax+1]
Definition:
runge_kutta_semiimplicit.icc:60
T
Expr1::float_type T
Definition:
field_expr.h:261