ROL
OrthogonalPolynomials.hpp
Go to the documentation of this file.
1#include<iostream>
2#include<iomanip>
3#include<cmath>
4#include"ROL_StdVector.hpp"
5#include"Teuchos_LAPACK.hpp"
6#include"LinearAlgebra.hpp"
7
8
9#ifndef __ORTHOGONAL_POLYNOMIALS__
10#define __ORTHOGONAL_POLYNOMIALS__
11
12
28template<class Real>
29void rec_jacobi( ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
30 const double alpha,
31 const double beta,
32 std::vector<Real> &a,
33 std::vector<Real> &b ) {
34
35 int N = a.size();
36 Real nu = (beta-alpha)/double(alpha+beta+2.0);
37 Real mu = pow(2.0,alpha+beta+1.0)*tgamma(alpha+1.0)*tgamma(beta+1.0)/tgamma(alpha+beta+2.0);
38 Real nab;
39 Real sqdif = pow(beta,2)-pow(alpha,2);
40
41 a[0] = nu;
42 b[0] = mu;
43
44 if(N>1){
45
46 for(int n=1;n<N;++n) {
47 nab = 2*n+alpha+beta;
48 a[n] = sqdif/(nab*(nab+2));
49 }
50
51 b[1] = 4.0*(alpha+1.0)*(beta+1.0)/(pow(alpha+beta+2.0,2)*(alpha+beta+3.0));
52
53 if(N>2) {
54 for(int n=2;n<N;++n) {
55 nab = 2*n+alpha+beta;
56 b[n] = 4.0*(n+alpha)*(n+beta)*n*(n+alpha+beta)/(nab*nab*(nab+1.0)*(nab-1.0));
57 }
58 }
59 }
60}
61
71template<class Real>
72void vandermonde( const std::vector<Real> &a,
73 const std::vector<Real> &b,
74 const std::vector<Real> &x,
75 std::vector<Real> &V) {
76 const int n = a.size();
77
78 for(int i=0;i<n;++i) {
79 // Zeroth polynomial is always 1
80 V[i] = 1.0;
81 // First polynomial is (x-a)
82 V[i+n] = x[i] - a[0];
83 }
84 for(int j=1;j<n-1;++j) {
85 for(int i=0;i<n;++i) {
86 V[i+(j+1)*n] = (x[i] - a[j])*V[i+j*n] - b[j]*V[i+(j-1)*n];
87 }
88 }
89}
90
91
103template<class Real>
104void gauss( ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
105 const std::vector<Real> &a,
106 const std::vector<Real> &b,
107 std::vector<Real> &x,
108 std::vector<Real> &w ) {
109 int INFO;
110 const int N = a.size();
111 const int LDZ = N;
112 const char COMPZ = 'I';
113
114 std::vector<Real> D(N,0.0);
115 std::vector<Real> E(N,0.0);
116 std::vector<Real> WORK(4*N,0.0);
117
118 // Column-stacked matrix of eigenvectors needed for weights
119 std::vector<Real> Z(N*N,0);
120
121 D = a;
122
123 for(int i=0;i<N-1;++i) {
124 E[i] = sqrt(b[i+1]);
125 }
126
127 lapack->STEQR(COMPZ,N,&D[0],&E[0],&Z[0],LDZ,&WORK[0],&INFO);
128
129 for(int i=0;i<N;++i){
130 x[i] = D[i];
131 w[i] = b[0]*pow(Z[N*i],2);
132 }
133
134}
135
136
137
150template<class Real>
151void rec_lobatto( ROL::Ptr<Teuchos::LAPACK<int,Real> > const lapack,
152 const double xl1,
153 const double xl2,
154 std::vector<Real> &a,
155 std::vector<Real> &b ) {
156 const int N = a.size()-1;
157
158 std::vector<Real> amod(N,0.0);
159 std::vector<Real> bmod(N-1,0.0);
160 std::vector<Real> en(N,0.0);
161 std::vector<Real> g(N,0.0);
162
163 // Nth canonical vector
164 en[N-1] = 1;
165
166 for(int i=0;i<N-1;++i) {
167 bmod[i] = sqrt(b[i+1]);
168 }
169
170 for(int i=0;i<N;++i) {
171 amod[i] = a[i]-xl1;
172 }
173
174 trisolve(lapack,bmod,amod,bmod,en,g);
175 Real g1 = g[N-1];
176
177 for(int i=0;i<N;++i) {
178 amod[i] = a[i]-xl2;
179 }
180
181 trisolve(lapack,bmod,amod,bmod,en,g);
182 Real g2 = g[N-1];
183
184 a[N] = (g1*xl2-g2*xl1)/(g1-g2);
185 b[N] = (xl2-xl1)/(g1-g2);
186
187}
188
189
190#endif
191
192
void trisolve(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &a, const std::vector< Real > &b, const std::vector< Real > &c, const std::vector< Real > &r, std::vector< Real > &x)
Solve a tridiagonal system.
void vandermonde(const std::vector< Real > &a, const std::vector< Real > &b, const std::vector< Real > &x, std::vector< Real > &V)
Construct the generalized Vandermonde matrix (in column stacked form) based upon the recurrence coeff...
void rec_jacobi(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const double alpha, const double beta, std::vector< Real > &a, std::vector< Real > &b)
Generate the Jacobi polynomial recursion coeffcients .
void gauss(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &a, const std::vector< Real > &b, std::vector< Real > &x, std::vector< Real > &w)
Compute the Gauss quadrature nodes and weights for the polynomials generated by the recurrence coeffi...
void rec_lobatto(ROL::Ptr< Teuchos::LAPACK< int, Real > > const lapack, const double xl1, const double xl2, std::vector< Real > &a, std::vector< Real > &b)
Modify the given recurrence coefficients so that the set of zeros of the maximal order polynomial inc...
Vector< Real > V