Nektar++
Polylib.cpp
Go to the documentation of this file.
1 #include "Polylib.h"
2 #include <cstdio>
3 #include <cmath>
4 #include <cfloat>
5 #include <complex>
6 
7 #include <boost/core/ignore_unused.hpp>
8 
9 #include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
10 /// Maximum number of iterations in polynomial defalation routine Jacobz
11 #define STOP 30
12 /// Precision tolerance for two points to be similar
13 #define EPS 100*DBL_EPSILON
14 /// return the sign(b)*a
15 #define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a))
16 
17 namespace Polylib {
18 
19  /// The following function is used to circumvent/reduce "Subtractive Cancellation"
20  /// The expression 1/dz is replaced by optinvsub(.,.)
21  /// Added on 26 April 2017
22  double optdiff(double xl, double xr)
23  {
24  double m_xln, m_xrn;
25  int m_expn;
26  int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
27 
28  if (fabs(xl-xr)<1.e-4){
29 
30  m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
31  m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
32  m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
33  m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
34  m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
35 
36  return powl(10.0L,-m_digits)*(m_xln-m_xrn);
37  }else{
38  return (xl-xr);
39  }
40  }
41 
42  double laginterp(double z, int j, const double *zj, int np)
43  {
44  double temp = 1.0;
45  for (int i=0; i<np; i++)
46  {
47  if (j != i)
48  {
49  temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
50  }
51  }
52  return temp;
53  }
54  /// Define whether to use polynomial deflation (1) or tridiagonal solver (0).
55 #define POLYNOMIAL_DEFLATION 0
56 
57 #ifdef POLYNOMIAL_DEFLATION
58  /// zero determination using Newton iteration with polynomial deflation
59 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
60 #else
61  /// zero determination using eigenvalues of tridiagaonl matrix
62 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
63 #endif
64 
65 
66 
67  /* local functions */
68  static void Jacobz (const int n, double *z, const double alpha,
69  const double beta);
70  // static void JacZeros (const int n, double *a, const double alpha,
71  // const double beta);
72  //static void TriQL (const int n, double *d, double *e);
73  static void TriQL(const int, double *,double *, double **);
74  double gammaF (const double);
75  double gammaFracGammaF(const int, const double, const int, const double);
76  static void RecCoeff(const int, double *, double *,const double,
77  const double);
78  void JKMatrix(int, double *, double *);
79  void chri1(int,double*,double*,double*,double*,double);
80 
81  /**
82  \brief Gauss-Jacobi zeros and weights.
83 
84  \li Generate \a np Gauss Jacobi zeros, \a z, and weights,\a w,
85  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
86  \f$,
87 
88  \li Exact for polynomials of order \a 2np-1 or less
89  */
90 
91  void zwgj (double *z, double *w, const int np, const double alpha,
92  const double beta)
93  {
94  int i;
95  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
96 
97  jacobz (np,z,alpha,beta);
98  jacobd (np,z,w,np,alpha,beta);
99 
100  fac = pow(two,apb + one) * gammaFracGammaF(np + 1, alpha, np + 1, 0.0)
101  * gammaFracGammaF(np + 1, beta , np + 1, apb);
102 
103  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
104 
105  return;
106  }
107 
108 
109  /**
110  \brief Gauss-Radau-Jacobi zeros and weights with end point at \a z=-1.
111 
112  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
113  associated with the polynomial \f$(1+z) P^{\alpha,\beta+1}_{np-1}(z)
114  \f$.
115 
116  \li Exact for polynomials of order \a 2np-2 or less
117  */
118 
119  void zwgrjm(double *z, double *w, const int np, const double alpha,
120  const double beta)
121  {
122 
123  if(np == 1){
124  z[0] = 0.0;
125  w[0] = 2.0;
126  }
127  else{
128  int i;
129  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
130 
131  z[0] = -one;
132  jacobz (np-1,z+1,alpha,beta+1);
133  jacobfd (np,z,w,NULL,np-1,alpha,beta);
134 
135  fac = pow(two,apb) * gammaFracGammaF(np, alpha, np, 0.0)
136  * gammaFracGammaF(np, beta , np + 1, apb);
137  fac /= (beta + np);
138 
139  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
140  w[0] *= (beta + one);
141  }
142 
143  return;
144  }
145 
146 
147  /**
148  \brief Gauss-Radau-Jacobi zeros and weights with end point at \a z=1
149 
150 
151  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
152  associated with the polynomial \f$(1-z) P^{\alpha+1,\beta}_{np-1}(z)
153  \f$.
154 
155  \li Exact for polynomials of order \a 2np-2 or less
156  */
157 
158  void zwgrjp(double *z, double *w, const int np, const double alpha,
159  const double beta)
160  {
161 
162  if(np == 1){
163  z[0] = 0.0;
164  w[0] = 2.0;
165  }
166  else{
167  int i;
168  double fac, one = 1.0, two = 2.0, apb = alpha + beta;
169 
170  jacobz (np-1,z,alpha+1,beta);
171  z[np-1] = one;
172  jacobfd (np,z,w,NULL,np-1,alpha,beta);
173 
174  fac = pow(two,apb) * gammaFracGammaF(np, alpha, np, 0.0)
175  * gammaFracGammaF(np, beta , np + 1, apb);
176  fac /= (alpha + np);
177 
178  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
179  w[np-1] *= (alpha + one);
180  }
181 
182  return;
183  }
184 
185 
186  /**
187  \brief Gauss-Lobatto-Jacobi zeros and weights with end point at \a z=-1,\a 1
188 
189 
190  \li Generate \a np Gauss-Lobatto-Jacobi points, \a z, and weights, \a w,
191  associated with polynomial \f$ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \f$
192  \li Exact for polynomials of order \a 2np-3 or less
193  */
194 
195  void zwglj(double *z, double *w, const int np, const double alpha,
196  const double beta)
197  {
198 
199  if( np == 1 ){
200  z[0] = 0.0;
201  w[0] = 2.0;
202  }
203  else if( np == 2 ){
204  z[0] = -1.0;
205  z[1] = 1.0;
206 
207  w[0] = 1.0;
208  w[1] = 1.0;
209  }
210  else{
211  int i;
212  double fac, one = 1.0, apb = alpha + beta, two = 2.0;
213 
214  z[0] = -one;
215  z[np-1] = one;
216  jacobz (np-2,z + 1,alpha + one,beta + one);
217  jacobfd (np,z,w,NULL,np-1,alpha,beta);
218 
219  fac = pow(two,apb + 1) * gammaFracGammaF(np, alpha, np, 0.0)
220  * gammaFracGammaF(np, beta , np + 1, apb);
221  fac /= (np-1);
222 
223  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
224  w[0] *= (beta + one);
225  w[np-1] *= (alpha + one);
226  }
227 
228  return;
229  }
230 
231  /**
232  \brief Gauss-Kronrod-Jacobi zeros and weights.
233 
234  \li Generate \a npt=2*np+1 Gauss-Kronrod Jacobi zeros, \a z, and weights,\a w,
235  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
236  \f$,
237 
238  \li Exact for polynomials of order \a 3np+1 or less
239  */
240  void zwgk(double *z, double *w, const int npt , const double alpha,
241  const double beta)
242  {
243 
244  int np = (npt-1)/2;
245 
246  int i,j;
247 
248  // number of kronrod points associated with the np gauss rule
249  int kpoints = 2*np + 1;
250 
251  // Define the number of required recurrence coefficents
252  int ncoeffs = (int)floor(3.0*(np+1)/2);
253 
254  // Define arrays for the recurrence coefficients
255  // We will use these arrays for the Kronrod results too, hence the
256  // reason for the size of the arrays
257  double *a = new double[kpoints];
258  double *b = new double[kpoints];
259 
260  // Initialize a and b to zero
261  for(i = 0; i < kpoints; i++)
262  {
263  a[i] = 0.0;
264  b[i] = 0.0;
265  }
266 
267  // Call the routine to calculate the recurrence coefficients
268  RecCoeff(ncoeffs,a,b,alpha,beta);
269 
270  // Call the routine to caluclate the jacobi-Kronrod matrix
271  JKMatrix(np,a,b);
272 
273  // Set up the identity matrix
274  double** zmatrix = new double*[kpoints];
275  for(i = 0; i < kpoints; i++)
276  {
277  zmatrix[i] = new double[kpoints];
278  for(j = 0; j < kpoints; j++)
279  {
280  zmatrix[i][j] = 0.0;
281  }
282  }
283  for(i = 0; i < kpoints; i++)
284  {
285  zmatrix[i][i] = 1.0;
286  }
287 
288  // Calculte the points and weights
289  TriQL(kpoints, a, b, zmatrix);
290 
291  for(i = 0; i < kpoints; i++)
292  {
293  z[i] = a[i];
294  w[i] = b[i];
295  }
296  delete[] a;
297  delete[] b;
298  for (i = 0; i < kpoints; i++)
299  {
300  delete[] zmatrix[i];
301  }
302  delete[] zmatrix;
303 
304  }
305 
306  /**
307  \brief Gauss-Radau-Kronrod-Jacobi zeros and weights.
308 
309  \li Generate \a npt=2*np Radau-Kronrod Jacobi zeros, \a z, and weights,\a w,
310  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
311  \f$,
312  */
313  void zwrk(double* z, double* w, const int npt ,const double alpha,
314  const double beta)
315  {
316 
317  int np = npt/2;
318 
319  if(np < 2)
320  {
321  fprintf(stderr,"too few points in formula\n");
322  return;
323  }
324 
325  double end0 = -1;
326 
327  int i,j;
328 
329  // number of kronrod points associated with the np gauss rule
330  int kpoints = 2*np;
331 
332  // Define the number of required recurrence coefficents
333  int ncoeffs = (int)ceil(3.0*np/2);
334 
335  // Define arrays for the recurrence coefficients
336  double *a = new double[ncoeffs+1];
337  double *b = new double[ncoeffs+1];
338 
339  // Initialize a and b to zero
340  for(i = 0; i < ncoeffs+1; i++)
341  {
342  a[i] = 0.0;
343  b[i] = 0.0;
344  }
345 
346  // Call the routine to calculate the recurrence coefficients
347  RecCoeff(ncoeffs,a,b,alpha,beta);
348 
349  double* a0 = new double[ncoeffs];
350  double* b0 = new double[ncoeffs];
351 
352  chri1(ncoeffs,a,b,a0,b0,end0);
353 
354  double s = b0[0]/fabs(b0[0]);
355  b0[0] = s*b0[0];
356 
357  // Finding the 2*np-1 gauss-kronrod points
358  double* z1 = new double[2*np-1];
359  double* w1 = new double[2*np-1];
360  for(i = 0; i < ncoeffs; i++)
361  {
362  z1[i] = a0[i];
363  w1[i] = b0[i];
364  }
365  JKMatrix(np-1,z1,w1);
366  // Set up the identity matrix
367  double** zmatrix = new double*[2*np-1];
368  for(i = 0; i < 2*np-1; i++)
369  {
370  zmatrix[i] = new double[2*np-1];
371  for(j = 0; j < 2*np-1; j++)
372  {
373  zmatrix[i][j] = 0.0;
374  }
375  }
376  for(i = 0; i < 2*np-1; i++)
377  {
378  zmatrix[i][i] = 1.0;
379  }
380 
381  // Calculate the points and weights
382  TriQL(2*np-1, z1, w1, zmatrix);
383 
384  double sumW1 = 0.0;
385  for(i = 0; i < 2*np-1; i++)
386  {
387  w1[i] = s*w1[i]/(z1[i]-end0);
388  sumW1 += w1[i];
389  }
390 
391  z[0] = end0;
392  w[0] = b[0]- sumW1;
393  for(i = 1; i < kpoints; i++)
394  {
395  z[i] = z1[i-1];
396  w[i] = w1[i-1];
397  }
398 
399 
400  delete[] a;
401  delete[] b;
402  delete[] a0;
403  delete[] b0;
404  delete[] z1;
405  delete[] w1;
406  for(i = 0; i < 2*np-1; i++)
407  {
408  delete[] zmatrix[i];
409  }
410  delete[] zmatrix;
411  }
412 
413  /**
414  \brief Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
415 
416  \li Generate \a npt=2*np-1 Lobatto-Kronrod Jacobi zeros, \a z, and weights,\a w,
417  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
418  \f$,
419  */
420  void zwlk(double* z, double* w, const int npt,
421  const double alpha, const double beta)
422  {
423 
424  int np = (npt+1)/2;
425 
426  if(np < 4)
427  {
428  fprintf (stderr,"too few points in formula\n");
429  return;
430  }
431 
432  double endl = -1;
433  double endr = 1;
434  int i,j;
435 
436  // number of kronrod points associated with the np gauss rule
437  int kpoints = 2*np-1;
438 
439  // Define the number of required recurrence coefficents
440  int ncoeffs = (int)ceil(3.0*np/2)-1;
441 
442  // Define arrays for the recurrence coefficients
443  double *a = new double[ncoeffs+1];
444  double *b = new double[ncoeffs+1];
445 
446  // Initialize a and b to zero
447  for(i = 0; i < ncoeffs+1; i++)
448  {
449  a[i] = 0.0;
450  b[i] = 0.0;
451  }
452 
453  // Call the routine to calculate the recurrence coefficients
454  RecCoeff(ncoeffs,a,b,alpha,beta);
455 
456 
457  double* a0 = new double[ncoeffs];
458  double* b0 = new double[ncoeffs];
459 
460  chri1(ncoeffs,a,b,a0,b0,endl);
461 
462  double* a1 = new double[ncoeffs-1];
463  double* b1 = new double[ncoeffs-1];
464 
465  chri1(ncoeffs-1,a0,b0,a1,b1,endr);
466 
467 
468  double s = b1[0]/fabs(b1[0]);
469  b1[0] = s*b1[0];
470 
471  // Finding the 2*np-1 gauss-kronrod points
472  double* z1 = new double[2*np-3];
473  double* w1 = new double[2*np-3];
474  for(i = 0; i < ncoeffs; i++)
475  {
476  z1[i] = a1[i];
477  w1[i] = b1[i];
478  }
479  JKMatrix(np-2,z1,w1);
480  // Set up the identity matrix
481  double** zmatrix = new double*[2*np-3];
482  for(i = 0; i < 2*np-3; i++)
483  {
484  zmatrix[i] = new double[2*np-3];
485  for(j = 0; j < 2*np-3; j++)
486  {
487  zmatrix[i][j] = 0.0;
488  }
489  }
490  for(i = 0; i < 2*np-3; i++)
491  {
492  zmatrix[i][i] = 1.0;
493  }
494 
495  // Calculate the points and weights
496  TriQL(2*np-3, z1, w1, zmatrix);
497 
498  double sumW1 = 0.0;
499  double sumW1Z1 = 0.0;
500  for(i = 0; i < 2*np-3; i++)
501  {
502  w1[i] = s*w1[i]/(z1[i]-endl)/(z1[i]-endr);
503  sumW1 += w1[i];
504  sumW1Z1 += z1[i]*w1[i];
505  }
506 
507  double c0 = b[0]-sumW1;
508  double c1 = a[0]*b[0]-sumW1Z1;
509 
510  z[0] = endl;
511  z[2*np-2] = endr;
512  w[0] = (c0*endr-c1)/(endr-endl);
513  w[2*np-2] = (c1-c0*endl)/(endr-endl);
514 
515  for(i = 1; i < kpoints-1; i++)
516  {
517  z[i] = z1[i-1];
518  w[i] = w1[i-1];
519  }
520  delete[] a;
521  delete[] b;
522  delete[] a0;
523  delete[] b0;
524  delete[] a1;
525  delete[] b1;
526  delete[] z1;
527  delete[] w1;
528  for(i = 0; i < 2*np-3; i++)
529  {
530  delete[] zmatrix[i];
531  }
532  delete[] zmatrix;
533  }
534 
535  /**
536  \brief Compute the Derivative Matrix and its transpose associated
537  with the Gauss-Jacobi zeros.
538 
539  \li Compute the derivative matrix, \a d, and its transpose, \a dt,
540  associated with the n_th order Lagrangian interpolants through the
541  \a np Gauss-Jacobi points \a z such that \n
542  \f$ \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
543 
544  */
545 
546  void Dgj(double *D, const double *z, const int np, const double alpha,
547  const double beta)
548  {
549 
550  double one = 1.0, two = 2.0;
551 
552  if (np <= 0){
553  D[0] = 0.0;
554  }
555  else{
556  int i,j;
557  double *pd;
558 
559  pd = (double *)malloc(np*sizeof(double));
560  jacobd(np,z,pd,np,alpha,beta);
561 
562  for (i = 0; i < np; i++){
563  for (j = 0; j < np; j++){
564 
565  if (i != j)
566  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
567  else
568  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
569  (two*(one - z[j]*z[j]));
570  }
571  }
572  free(pd);
573  }
574  return;
575  }
576 
577 
578  /**
579  \brief Compute the Derivative Matrix and its transpose associated
580  with the Gauss-Radau-Jacobi zeros with a zero at \a z=-1.
581 
582  \li Compute the derivative matrix, \a d, associated with the n_th
583  order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
584  points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
585  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
586 
587  */
588 
589  void Dgrjm(double *D, const double *z, const int np, const double alpha,
590  const double beta)
591  {
592 
593  if (np <= 0){
594  D[0] = 0.0;
595  }
596  else{
597  int i, j;
598  double one = 1.0, two = 2.0;
599  double *pd;
600 
601  pd = (double *)malloc(np*sizeof(double));
602 
603  pd[0] = pow(-one,np-1) * gammaFracGammaF(np + 1, beta, np, 0.0);
604  pd[0] /= gammaF(beta+two);
605  jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
606  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
607 
608  for (i = 0; i < np; i++) {
609  for (j = 0; j < np; j++){
610  if (i != j)
611  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
612  else {
613  if(j == 0)
614  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
615  (two*(beta + two));
616  else
617  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
618  (two*(one - z[j]*z[j]));
619  }
620  }
621  }
622  free(pd);
623  }
624 
625  return;
626  }
627 
628 
629  /**
630  \brief Compute the Derivative Matrix associated with the
631  Gauss-Radau-Jacobi zeros with a zero at \a z=1.
632 
633  \li Compute the derivative matrix, \a d, associated with the n_th
634  order Lagrangian interpolants through the \a np Gauss-Radau-Jacobi
635  points \a z such that \n \f$ \frac{du}{dz}(z[i]) =
636  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
637  */
638 
639  void Dgrjp(double *D, const double *z, const int np, const double alpha,
640  const double beta)
641  {
642 
643  if (np <= 0){
644  D[0] = 0.0;
645  }
646  else{
647  int i, j;
648  double one = 1.0, two = 2.0;
649  double *pd;
650 
651  pd = (double *)malloc(np*sizeof(double));
652 
653 
654  jacobd(np-1,z,pd,np-1,alpha+1,beta);
655  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
656  pd[np-1] = -gammaFracGammaF(np + 1, alpha, np, 0.0);
657  pd[np-1] /= gammaF(alpha+two);
658 
659  for (i = 0; i < np; i++) {
660  for (j = 0; j < np; j++){
661  if (i != j)
662  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
663  else {
664  if(j == np-1)
665  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
666  (two*(alpha + two));
667  else
668  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
669  (two*(one - z[j]*z[j]));
670  }
671  }
672  }
673  free(pd);
674  }
675 
676  return;
677  }
678 
679  /**
680  \brief Compute the Derivative Matrix associated with the
681  Gauss-Lobatto-Jacobi zeros.
682 
683  \li Compute the derivative matrix, \a d, associated with the n_th
684  order Lagrange interpolants through the \a np
685  Gauss-Lobatto-Jacobi points \a z such that \n \f$
686  \frac{du}{dz}(z[i]) = \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$
687 
688  */
689 
690  void Dglj(double *D, const double *z, const int np, const double alpha,
691  const double beta)
692  {
693  if (np <= 1){
694  D[0] = 0.0;
695  }
696  else{
697  int i, j;
698  double one = 1.0, two = 2.0;
699  double *pd;
700 
701  pd = (double *)malloc(np*sizeof(double));
702 
703  pd[0] = two*pow(-one,np) * gammaFracGammaF(np, beta, np - 1, 0.0);
704  pd[0] /= gammaF(beta + two);
705  jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
706  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
707  pd[np-1] = -two * gammaFracGammaF(np, alpha, np - 1, 0.0);
708  pd[np-1] /= gammaF(alpha + two);
709 
710  for (i = 0; i < np; i++) {
711  for (j = 0; j < np; j++){
712  if (i != j)
713  D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i]));
714  else {
715  if (j == 0)
716  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
717  else if (j == np-1)
718  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
719  else
720  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
721  (two*(one - z[j]*z[j]));
722  }
723  }
724  }
725  free(pd);
726  }
727 
728  return;
729  }
730 
731 
732  /**
733  \brief Compute the value of the \a i th Lagrangian interpolant through
734  the \a np Gauss-Jacobi points \a zgj at the arbitrary location \a z.
735 
736  \li \f$ -1 \leq z \leq 1 \f$
737 
738  \li Uses the defintion of the Lagrangian interpolant:\n
739  %
740  \f$ \begin{array}{rcl}
741  h_j(z) = \left\{ \begin{array}{ll}
742  \displaystyle \frac{P_{np}^{\alpha,\beta}(z)}
743  {[P_{np}^{\alpha,\beta}(z_j)]^\prime
744  (z-z_j)} & \mbox{if $z \ne z_j$}\\
745  & \\
746  1 & \mbox{if $z=z_j$}
747  \end{array}
748  \right.
749  \end{array} \f$
750  */
751 
752  double hgj (const int i, const double z, const double *zgj,
753  const int np, const double alpha, const double beta)
754  {
755  boost::ignore_unused(alpha, beta);
756  double zi, dz;
757 
758  zi = *(zgj+i);
759  dz = z-zi;
760  if (fabs(dz) < EPS) return 1.0;
761 
762  return laginterp(z, i, zgj, np);
763 
764  }
765 
766  /**
767  \brief Compute the value of the \a i th Lagrangian interpolant through the
768  \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
769  \a z. This routine assumes \a zgrj includes the point \a -1.
770 
771  \li \f$ -1 \leq z \leq 1 \f$
772 
773  \li Uses the defintion of the Lagrangian interpolant:\n
774  %
775  \f$ \begin{array}{rcl}
776  h_j(z) = \left\{ \begin{array}{ll}
777  \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)}
778  {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime +
779  P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
780  & \\
781  1 & \mbox{if $z=z_j$}
782  \end{array}
783  \right.
784  \end{array} \f$
785  */
786 
787  double hgrjm (const int i, const double z, const double *zgrj, const int np,
788  const double alpha, const double beta)
789  {
790  boost::ignore_unused(alpha, beta);
791 
792  double zi, dz;
793 
794  zi = *(zgrj+i);
795  dz = z-zi;
796  if (fabs(dz) < EPS) return 1.0;
797 
798  return laginterp(z, i, zgrj, np);
799  }
800 
801 
802  /**
803  \brief Compute the value of the \a i th Lagrangian interpolant through the
804  \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
805  \a z. This routine assumes \a zgrj includes the point \a +1.
806 
807  \li \f$ -1 \leq z \leq 1 \f$
808 
809  \li Uses the defintion of the Lagrangian interpolant:\n
810  %
811  \f$ \begin{array}{rcl}
812  h_j(z) = \left\{ \begin{array}{ll}
813  \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)}
814  {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime -
815  P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
816  & \\
817  1 & \mbox{if $z=z_j$}
818  \end{array}
819  \right.
820  \end{array} \f$
821  */
822 
823  double hgrjp (const int i, const double z, const double *zgrj, const int np,
824  const double alpha, const double beta)
825  {
826  boost::ignore_unused(alpha, beta);
827  double zi, dz;
828 
829  zi = *(zgrj+i);
830  dz = z-zi;
831  if (fabs(dz) < EPS) return 1.0;
832 
833  return laginterp(z, i, zgrj, np);
834  }
835 
836 
837  /**
838  \brief Compute the value of the \a i th Lagrangian interpolant through the
839  \a np Gauss-Lobatto-Jacobi points \a zgrj at the arbitrary location
840  \a z.
841 
842  \li \f$ -1 \leq z \leq 1 \f$
843 
844  \li Uses the defintion of the Lagrangian interpolant:\n
845  %
846  \f$ \begin{array}{rcl}
847  h_j(z) = \left\{ \begin{array}{ll}
848  \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)}
849  {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime -
850  2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\
851  & \\
852  1 & \mbox{if $z=z_j$}
853  \end{array}
854  \right.
855  \end{array} \f$
856  */
857 
858  double hglj (const int i, const double z, const double *zglj, const int np,
859  const double alpha, const double beta)
860  {
861  boost::ignore_unused(alpha, beta);
862 
863  double zi, dz;
864 
865  zi = *(zglj+i);
866  dz = z-zi;
867  if (fabs(dz) < EPS) return 1.0;
868 
869  return laginterp(z, i, zglj, np);
870 
871  }
872 
873 
874  /**
875  \brief Interpolation Operator from Gauss-Jacobi points to an
876  arbitrary distribution at points \a zm
877 
878  \li Computes the one-dimensional interpolation matrix, \a im, to
879  interpolate a function from at Gauss-Jacobi distribution of \a nz
880  zeros \a zgrj to an arbitrary distribution of \a mz points \a zm, i.e.\n
881  \f$
882  u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j])
883  \f$
884 
885  */
886 
887  void Imgj(double *im, const double *zgj, const double *zm, const int nz,
888  const int mz,const double alpha, const double beta){
889  double zp;
890  int i, j;
891 
892  for (i = 0; i < nz; ++i) {
893  for (j = 0; j < mz; ++j)
894  {
895  zp = zm[j];
896  im [i*mz+j] = hgj(i, zp, zgj, nz, alpha, beta);
897  }
898  }
899 
900  return;
901  }
902 
903  /**
904  \brief Interpolation Operator from Gauss-Radau-Jacobi points
905  (including \a z=-1) to an arbitrary distrubtion at points \a zm
906 
907  \li Computes the one-dimensional interpolation matrix, \a im, to
908  interpolate a function from at Gauss-Radau-Jacobi distribution of
909  \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
910  distribution of \a mz points \a zm, i.e.
911  \n
912  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
913 
914  */
915 
916  void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz,
917  const int mz, const double alpha, const double beta)
918  {
919  double zp;
920  int i, j;
921 
922  for (i = 0; i < nz; i++) {
923  for (j = 0; j < mz; j++)
924  {
925  zp = zm[j];
926  im [i*mz+j] = hgrjm(i, zp, zgrj, nz, alpha, beta);
927  }
928  }
929 
930  return;
931  }
932 
933  /**
934  \brief Interpolation Operator from Gauss-Radau-Jacobi points
935  (including \a z=1) to an arbitrary distrubtion at points \a zm
936 
937  \li Computes the one-dimensional interpolation matrix, \a im, to
938  interpolate a function from at Gauss-Radau-Jacobi distribution of
939  \a nz zeros \a zgrj (where \a zgrj[nz-1]=1) to an arbitrary
940  distribution of \a mz points \a zm, i.e.
941  \n
942  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
943 
944  */
945 
946  void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz,
947  const int mz,const double alpha, const double beta)
948  {
949  double zp;
950  int i, j;
951 
952  for (i = 0; i < nz; i++) {
953  for (j = 0; j < mz; j++)
954  {
955  zp = zm[j];
956  im [i*mz+j] = hgrjp(i, zp, zgrj, nz, alpha, beta);
957  }
958  }
959 
960  return;
961  }
962 
963 
964  /**
965  \brief Interpolation Operator from Gauss-Lobatto-Jacobi points
966  to an arbitrary distrubtion at points \a zm
967 
968  \li Computes the one-dimensional interpolation matrix, \a im, to
969  interpolate a function from at Gauss-Lobatto-Jacobi distribution of
970  \a nz zeros \a zgrj (where \a zgrj[0]=-1) to an arbitrary
971  distribution of \a mz points \a zm, i.e.
972  \n
973  \f$ u(zm[i]) = \sum_{j=0}^{nz-1} im[i*nz+j] \ u(zgj[j]) \f$
974 
975  */
976 
977  void Imglj(double *im, const double *zglj, const double *zm, const int nz,
978  const int mz, const double alpha, const double beta)
979  {
980  double zp;
981  int i, j;
982 
983  for (i = 0; i < nz; i++) {
984  for (j = 0; j < mz; j++)
985  {
986  zp = zm[j];
987  im[i*mz+j] = hglj(i, zp, zglj, nz, alpha, beta);
988  }
989  }
990 
991  return;
992  }
993 
994  /**
995  \brief Routine to calculate Jacobi polynomials, \f$
996  P^{\alpha,\beta}_n(z) \f$, and their first derivative, \f$
997  \frac{d}{dz} P^{\alpha,\beta}_n(z) \f$.
998 
999  \li This function returns the vectors \a poly_in and \a poly_d
1000  containing the value of the \f$ n^th \f$ order Jacobi polynomial
1001  \f$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \f$ and its
1002  derivative at the \a np points in \a z[i]
1003 
1004  - If \a poly_in = NULL then only calculate derivatice
1005 
1006  - If \a polyd = NULL then only calculate polynomial
1007 
1008  - To calculate the polynomial this routine uses the recursion
1009  relationship (see appendix A ref [4]) :
1010  \f$ \begin{array}{rcl}
1011  P^{\alpha,\beta}_0(z) &=& 1 \\
1012  P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\
1013  a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z)
1014  P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\
1015  a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\
1016  a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2) \\
1017  a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1)
1018  (2n + \alpha + \beta + 2) \\
1019  a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2)
1020  \end{array} \f$
1021 
1022  - To calculate the derivative of the polynomial this routine uses
1023  the relationship (see appendix A ref [4]) :
1024  \f$ \begin{array}{rcl}
1025  b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z)
1026  + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\
1027  b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\
1028  b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\
1029  b^3_n(z) &=& 2(n+\alpha)(n+\beta)
1030  \end{array} \f$
1031 
1032  - Note the derivative from this routine is only valid for -1 < \a z < 1.
1033  */
1034  void jacobfd(const int np, const double *z, double *poly_in, double *polyd,
1035  const int n, const double alpha, const double beta){
1036  int i;
1037  double zero = 0.0, one = 1.0, two = 2.0;
1038 
1039  if(!np)
1040  return;
1041 
1042  if(n == 0){
1043  if(poly_in)
1044  for(i = 0; i < np; ++i)
1045  poly_in[i] = one;
1046  if(polyd)
1047  for(i = 0; i < np; ++i)
1048  polyd[i] = zero;
1049  }
1050  else if (n == 1){
1051  if(poly_in)
1052  for(i = 0; i < np; ++i)
1053  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1054  if(polyd)
1055  for(i = 0; i < np; ++i)
1056  polyd[i] = 0.5*(alpha + beta + two);
1057  }
1058  else{
1059  int k;
1060  double a1,a2,a3,a4;
1061  double two = 2.0, apb = alpha + beta;
1062  double *poly, *polyn1,*polyn2;
1063 
1064  if(poly_in){ // switch for case of no poynomial function return
1065  polyn1 = (double *)malloc(2*np*sizeof(double));
1066  polyn2 = polyn1+np;
1067  poly = poly_in;
1068  }
1069  else{
1070  polyn1 = (double *)malloc(3*np*sizeof(double));
1071  polyn2 = polyn1+np;
1072  poly = polyn2+np;
1073  }
1074 
1075  for(i = 0; i < np; ++i){
1076  polyn2[i] = one;
1077  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
1078  }
1079 
1080  for(k = 2; k <= n; ++k){
1081  a1 = two*k*(k + apb)*(two*k + apb - two);
1082  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
1083  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
1084  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
1085 
1086  a2 /= a1;
1087  a3 /= a1;
1088  a4 /= a1;
1089 
1090  for(i = 0; i < np; ++i){
1091  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
1092  polyn2[i] = polyn1[i];
1093  polyn1[i] = poly [i];
1094  }
1095  }
1096 
1097  if(polyd){
1098  a1 = n*(alpha - beta);
1099  a2 = n*(two*n + alpha + beta);
1100  a3 = two*(n + alpha)*(n + beta);
1101  a4 = (two*n + alpha + beta);
1102  a1 /= a4; a2 /= a4; a3 /= a4;
1103 
1104  // note polyn2 points to polyn1 at end of poly iterations
1105  for(i = 0; i < np; ++i){
1106  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
1107  polyd[i] /= (one - z[i]*z[i]);
1108  }
1109  }
1110 
1111  free(polyn1);
1112  }
1113 
1114  return;
1115  }
1116 
1117 
1118  /**
1119  \brief Calculate the derivative of Jacobi polynomials
1120 
1121  \li Generates a vector \a poly of values of the derivative of the
1122  \a n th order Jacobi polynomial \f$ P^(\alpha,\beta)_n(z)\f$ at the
1123  \a np points \a z.
1124 
1125  \li To do this we have used the relation
1126  \n
1127  \f$ \frac{d}{dz} P^{\alpha,\beta}_n(z)
1128  = \frac{1}{2} (\alpha + \beta + n + 1) P^{\alpha,\beta}_n(z) \f$
1129 
1130  \li This formulation is valid for \f$ -1 \leq z \leq 1 \f$
1131 
1132  */
1133 
1134  void jacobd(const int np, const double *z, double *polyd, const int n,
1135  const double alpha, const double beta)
1136  {
1137  int i;
1138  double one = 1.0;
1139  if(n == 0)
1140  for(i = 0; i < np; ++i) polyd[i] = 0.0;
1141  else{
1142  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
1143  jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
1144  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (double)n + one);
1145  }
1146  return;
1147  }
1148 
1149 
1150  /**
1151  \brief Calculate the Gamma function , \f$ \Gamma(n)\f$, for integer
1152  values and halves.
1153 
1154  Determine the value of \f$\Gamma(n)\f$ using:
1155 
1156  \f$ \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$
1157 
1158  where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
1159  */
1160 
1161  double gammaF(const double x){
1162  double gamma = 1.0;
1163 
1164  if (x == -0.5) gamma = -2.0*sqrt(M_PI);
1165  else if (!x) return gamma;
1166  else if ((x-(int)x) == 0.5){
1167  int n = (int) x;
1168  double tmp = x;
1169 
1170  gamma = sqrt(M_PI);
1171  while(n--){
1172  tmp -= 1.0;
1173  gamma *= tmp;
1174  }
1175  }
1176  else if ((x-(int)x) == 0.0){
1177  int n = (int) x;
1178  double tmp = x;
1179 
1180  while(--n){
1181  tmp -= 1.0;
1182  gamma *= tmp;
1183  }
1184  }
1185  else
1186  fprintf(stderr,"%lf is not of integer or half order\n",x);
1187  return gamma;
1188  }
1189 
1190  /**
1191  \brief Calculate fraction of two Gamma functions, \f$ \Gamma(x+\alpha)/\Gamma(y+\beta) \f$,
1192  for integer values and halves.
1193 
1194  Determine the value of \f$\Gamma(n)\f$ using:
1195 
1196  \f$ \Gamma(n) = (n-1)! \mbox{ or } \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$
1197 
1198  where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
1199 
1200  Attempts simplification in cases like \f$ \Gamma(x+1)/\Gamma(x-1) = x*(x-1) \f$
1201  */
1202 
1203  double gammaFracGammaF(const int x, const double alpha, const int y, const double beta){
1204  double gamma = 1.0;
1205  double halfa = fabs(alpha - int(alpha));
1206  double halfb = fabs(beta - int(beta ));
1207  if (halfa == 0.0 && halfb == 0.0)// integer value
1208  {
1209  int X = x+alpha;
1210  int Y = y+beta;
1211  if (X > Y)
1212  {
1213  for(int tmp=X-1; tmp > Y-1; tmp-=1)
1214  gamma *= tmp;
1215  }
1216  else if (Y > X)
1217  {
1218  for(int tmp=Y-1; tmp > X-1; tmp-=1)
1219  gamma *= tmp;
1220  gamma = 1. / gamma;
1221  }
1222  }
1223  else if (halfa == 0.5 && halfb == 0.5)// both are halves
1224  {
1225  double X = x + alpha;
1226  double Y = y + beta;
1227  if( X > Y )
1228  {
1229  for(int tmp=int(X); tmp > int(Y); tmp-=1)
1230  gamma *= tmp-0.5;
1231  }
1232  else if( Y > X)
1233  {
1234  for(int tmp=int(Y); tmp > int(X); tmp-=1)
1235  gamma *= tmp-0.5;
1236  gamma = 1. / gamma;
1237  }
1238  }
1239  else
1240  {
1241  double X = x + alpha;
1242  double Y = y + beta;
1243  while (X>1 || Y>1)
1244  {
1245  if (X>1)
1246  {
1247  gamma *= X-1.;
1248  X -= 1.;
1249  }
1250  if (Y>1)
1251  {
1252  gamma /= Y-1.;
1253  Y -= 1.;
1254  }
1255  }
1256  if (X==0.5)
1257  {
1258  gamma *= sqrt(M_PI);
1259  }
1260  else if (Y==0.5)
1261  {
1262  gamma /= sqrt(M_PI);
1263  }
1264  else
1265  {
1266  fprintf(stderr,"%lf or %lf is not of integer or half order\n",X, Y);
1267  }
1268  }
1269 
1270  return gamma;
1271  }
1272 
1273  /**
1274  \brief Calculate the \a n zeros, \a z, of the Jacobi polynomial, i.e.
1275  \f$ P_n^{\alpha,\beta}(z) = 0 \f$
1276 
1277  This routine is only value for \f$( \alpha > -1, \beta > -1)\f$
1278  and uses polynomial deflation in a Newton iteration
1279  */
1280 
1281  static void Jacobz(const int n, double *z, const double alpha,
1282  const double beta){
1283  int i,j,k;
1284  double dth = M_PI/(2.0*(double)n);
1285  double poly,pder,rlast=0.0;
1286  double sum,delr,r;
1287  double one = 1.0, two = 2.0;
1288 
1289  if(!n)
1290  return;
1291 
1292  for(k = 0; k < n; ++k){
1293  r = -cos((two*(double)k + one) * dth);
1294  if(k) r = 0.5*(r + rlast);
1295 
1296  for(j = 1; j < STOP; ++j){
1297  jacobfd(1,&r,&poly, &pder, n, alpha, beta);
1298 
1299  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
1300 
1301  delr = -poly / (pder - sum * poly);
1302  r += delr;
1303  if( fabs(delr) < EPS ) break;
1304  }
1305  z[k] = r;
1306  rlast = r;
1307  }
1308  return;
1309  }
1310 
1311 
1312  /**
1313  \brief Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal
1314  matrix from the three term recurrence relationship.
1315 
1316  Set up a symmetric tridiagonal matrix
1317 
1318  \f$ \left [ \begin{array}{ccccc}
1319  a[0] & b[0] & & & \\
1320  b[0] & a[1] & b[1] & & \\
1321  0 & \ddots & \ddots & \ddots & \\
1322  & & \ddots & \ddots & b[n-2] \\
1323  & & & b[n-2] & a[n-1] \end{array} \right ] \f$
1324 
1325  Where the coefficients a[n], b[n] come from the recurrence relation
1326 
1327  \f$ b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1} p_{j-2}(z) \f$
1328 
1329  where \f$ j=n+1\f$ and \f$p_j(z)\f$ are the Jacobi (normalized)
1330  orthogonal polynomials \f$ \alpha,\beta > -1\f$( integer values and
1331  halves). Since the polynomials are orthonormalized, the tridiagonal
1332  matrix is guaranteed to be symmetric. The eigenvalues of this
1333  matrix are the zeros of the Jacobi polynomial.
1334  */
1335 
1336  void JacZeros(const int n, double *a, double*b, const double alpha,
1337  const double beta){
1338 
1339  int i,j;
1340  RecCoeff(n,a,b,alpha,beta);
1341 
1342  double **z = new double*[n];
1343  for(i = 0; i < n; i++)
1344  {
1345  z[i] = new double[n];
1346  for(j = 0; j < n; j++)
1347  {
1348  z[i][j] = 0.0;
1349  }
1350  }
1351  for(i = 0; i < n; i++)
1352  {
1353  z[i][i] = 1.0;
1354  }
1355 
1356  // find eigenvalues and eigenvectors
1357  TriQL(n, a, b,z);
1358 
1359  delete[] z;
1360  return;
1361  }
1362 
1363  /**
1364  \brief The routine finds the recurrence coefficients \a a and
1365  \a b of the orthogonal polynomials
1366  */
1367  static void RecCoeff(const int n, double *a, double *b,const double alpha,
1368  const double beta){
1369 
1370  int i;
1371  double apb, apbi,a2b2;
1372 
1373  if(!n)
1374  return;
1375 
1376  // generate normalised terms
1377  apb = alpha + beta;
1378  apbi = 2.0 + apb;
1379 
1380  b[0] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi); //MuZero
1381  a[0] = (beta-alpha)/apbi;
1382  b[1] = (4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
1383 
1384  a2b2 = beta*beta-alpha*alpha;
1385 
1386  for(i = 1; i < n-1; i++){
1387  apbi = 2.0*(i+1) + apb;
1388  a[i] = a2b2/((apbi-2.0)*apbi);
1389  b[i+1] = (4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
1390  ((apbi*apbi-1)*apbi*apbi));
1391  }
1392 
1393  apbi = 2.0*n + apb;
1394  a[n-1] = a2b2/((apbi-2.0)*apbi);
1395 
1396  }
1397 
1398 
1399  /** \brief QL algorithm for symmetric tridiagonal matrix
1400 
1401  This subroutine is a translation of an algol procedure,
1402  num. math. \b 12, 377-383(1968) by martin and wilkinson, as modified
1403  in num. math. \b 15, 450(1970) by dubrulle. Handbook for
1404  auto. comp., vol.ii-linear algebra, 241-248(1971). This is a
1405  modified version from numerical recipes.
1406 
1407  This subroutine finds the eigenvalues and first components of the
1408  eigenvectors of a symmetric tridiagonal matrix by the implicit QL
1409  method.
1410 
1411  on input:
1412  - n is the order of the matrix;
1413  - d contains the diagonal elements of the input matrix;
1414  - e contains the subdiagonal elements of the input matrix
1415  in its first n-2 positions.
1416  - z is the n by n identity matrix
1417 
1418  on output:
1419 
1420  - d contains the eigenvalues in ascending order.
1421  - e contains the weight values - modifications of the first component
1422  of normalised eigenvectors
1423  */
1424 
1425  static void TriQL(const int n, double *d,double *e, double **z){
1426  int m,l,iter,i,k;
1427  double s,r,p,g,f,dd,c,b;
1428 
1429  double MuZero = e[0];
1430 
1431  // Renumber the elements of e
1432  for(i = 0; i < n-1; i++)
1433  {
1434  e[i] = sqrt(e[i+1]);
1435  }
1436  e[n-1] = 0.0;
1437 
1438 
1439  for (l=0;l<n;l++) {
1440  iter=0;
1441  do {
1442  for (m=l;m<n-1;m++) {
1443  dd=fabs(d[m])+fabs(d[m+1]);
1444  if (fabs(e[m])+dd == dd) break;
1445  }
1446  if (m != l) {
1447  if (iter++ == STOP){
1448  fprintf(stderr,"triQL: Too many iterations in TQLI");
1449  exit(1);
1450  }
1451  g=(d[l+1]-d[l])/(2.0*e[l]);
1452  r=sqrt((g*g)+1.0);
1453  g=d[m]-d[l]+e[l]/(g+sign(r,g));
1454  s=c=1.0;
1455  p=0.0;
1456  for (i=m-1;i>=l;i--) {
1457  f=s*e[i];
1458  b=c*e[i];
1459  if (fabs(f) >= fabs(g)) {
1460  c=g/f;
1461  r=sqrt((c*c)+1.0);
1462  e[i+1]=f*r;
1463  c *= (s=1.0/r);
1464  } else {
1465  s=f/g;
1466  r=sqrt((s*s)+1.0);
1467  e[i+1]=g*r;
1468  s *= (c=1.0/r);
1469  }
1470  g=d[i+1]-p;
1471  r=(d[i]-g)*s+2.0*c*b;
1472  p=s*r;
1473  d[i+1]=g+p;
1474  g=c*r-b;
1475 
1476  // Calculate the eigenvectors
1477  for(k = 0; k < n; k++)
1478  {
1479  f = z[k][i+1];
1480  z[k][i+1] = s*z[k][i] + c*f;
1481  z[k][i] = c*z[k][i] - s*f;
1482  }
1483 
1484  }
1485  d[l]=d[l]-p;
1486  e[l]=g;
1487  e[m]=0.0;
1488  }
1489  } while (m != l);
1490  }
1491 
1492  // order eigenvalues
1493  // Since we only need the first component of the eigenvectors
1494  // to calcualte the weight, we only swap the first components
1495  for(i = 0; i < n-1; ++i){
1496  k = i;
1497  p = d[i];
1498  for(l = i+1; l < n; ++l)
1499  if (d[l] < p) {
1500  k = l;
1501  p = d[l];
1502  }
1503  d[k] = d[i];
1504  d[i] = p;
1505 
1506  double temp = z[0][k];
1507  z[0][k] = z[0][i];
1508  z[0][i] = temp;
1509  }
1510 
1511  // Calculate the weights
1512  for(i =0 ; i < n; i++)
1513  {
1514  e[i] = MuZero*z[0][i]*z[0][i];
1515  }
1516  }
1517 
1518  /**
1519  \brief Calcualtes the Jacobi-kronrod matrix by determining the
1520  \a a and \b coefficients.
1521 
1522  The first \a 3n+1 coefficients are already known
1523 
1524  For more information refer to:
1525  "Dirk P. Laurie, Calcualtion of Gauss-Kronrod quadrature rules"
1526  */
1527  void JKMatrix(int n, double *a, double *b)
1528  {
1529  int i,j,k,m;
1530  // Working storage
1531  int size = (int)floor(n/2.0)+2;
1532  double *s = new double[size];
1533  double *t = new double[size];
1534 
1535  // Initialize s and t to zero
1536  for(i = 0; i < size; i++)
1537  {
1538  s[i] = 0.0;
1539  t[i] = 0.0;
1540  }
1541 
1542  t[1] = b[n+1];
1543  for(m = 0; m <= n-2; m++)
1544  {
1545  double u = 0.0;
1546  for(k = (int)floor((m+1)/2.0); k >= 0; k--)
1547  {
1548  int l = m-k;
1549  u = u+(a[k+n+1]-a[l])*t[k+1] + b[k+n+1]*s[k] - b[l]*s[k+1];
1550  s[k+1] = u;
1551  }
1552 
1553  // Swap the contents of s and t
1554  double *hold = s;
1555  s = t;
1556  t = hold;
1557  }
1558 
1559 
1560  for(j = (int)floor(n/2.0); j >= 0; j--)
1561  {
1562  s[j+1] = s[j];
1563  }
1564 
1565  for(m = n-1; m <= 2*n-3; m++)
1566  {
1567  double u = 0;
1568  for(k = m+1-n; k <= floor((m-1)/2.0); k++)
1569  {
1570  int l = m-k;
1571  j = n-1-l;
1572  u = u-(a[k+n+1]-a[l])*t[j+1] - b[k+n+1]*s[j+1] + b[l]*s[j+2];
1573  s[j+1] = u;
1574  }
1575 
1576  if(m%2 == 0)
1577  {
1578  k = m/2;
1579  a[k+n+1] = a[k] + (s[j+1]-b[k+n+1]*s[j+2])/t[j+2];
1580 
1581  }else
1582  {
1583  k = (m+1)/2;
1584  b[k+n+1] = s[j+1]/s[j+2];
1585  }
1586 
1587 
1588  // Swap the contents of s and t
1589  double *hold = s;
1590  s = t;
1591  t = hold;
1592  }
1593 
1594  a[2*n ] = a[n-1]-b[2*n]*s[1]/t[1];
1595 
1596  }
1597 
1598  /**
1599  \brief
1600  Given a weight function \f$w(t)\f$ through the first \a n+1
1601  coefficients \a a and \a b of its orthogonal polynomials
1602  this routine generates the first \a n recurrence coefficients for the orthogonal
1603  polynomials relative to the modified weight function \f$(t-z)w(t)\f$.
1604 
1605  The result will be placed in the array \a a0 and \a b0.
1606  */
1607 
1608  void chri1(int n, double* a, double* b, double* a0,
1609  double* b0,double z)
1610  {
1611 
1612  double q = ceil(3.0*n/2);
1613  int size = (int)q+1;
1614  if(size < n+1)
1615  {
1616  fprintf(stderr,"input arrays a and b are too short\n");
1617  }
1618  double* r = new double[n+1];
1619  r[0] = z - a[0];
1620  r[1] = z - a[1] - b[1]/r[0];
1621  a0[0] = a[1] + r[1] - r[0];
1622  b0[0] = -r[0]*b[0];
1623 
1624  if(n == 1)
1625  {
1626  delete[] r;
1627  return;
1628  }
1629  int k = 0;
1630  for(k = 1; k < n; k++)
1631  {
1632  r[k+1] = z - a[k+1] - b[k+1]/r[k];
1633  a0[k] = a[k+1] + r[k+1] - r[k];
1634  b0[k] = b[k] * r[k]/r[k-1];
1635  }
1636  delete[] r;
1637 
1638  }
1639 
1640  /**
1641 
1642  \brief
1643 
1644  Calcualte the bessel function of the first kind with complex double input y.
1645  Taken from Numerical Recipies in C
1646 
1647  Returns a complex double
1648  */
1649 
1650 
1651  std::complex<Nektar::NekDouble> ImagBesselComp(int n,std::complex<Nektar::NekDouble> y)
1652  {
1653  std::complex<Nektar::NekDouble> z (1.0,0.0);
1654  std::complex<Nektar::NekDouble> zbes (1.0,0.0);
1655  std::complex<Nektar::NekDouble> zarg;
1656  Nektar::NekDouble tol = 1e-15;
1657  int maxit = 10000;
1658  int i = 1;
1659 
1660  zarg = -0.25*y*y;
1661 
1662  while (abs(z) > tol && i <= maxit){
1663  z = z*(1.0/i/(i+n)*zarg);
1664  if (abs(z) <= tol) break;
1665  zbes = zbes + z;
1666  i++;
1667  }
1668  zarg = 0.5*y;
1669  for (i=1;i<=n;i++){
1670  zbes = zbes*zarg;
1671  }
1672  return zbes;
1673 
1674  }
1675 } // end of namespace
1676 
jacobz
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
Definition: Polylib.cpp:59
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
EPS
#define EPS
Precision tolerance for two points to be similar.
Definition: Polylib.cpp:13
STOP
#define STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
Definition: Polylib.cpp:11
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Polylib
The namespace associated with the the Polylib library (Polylib introduction)
Definition: Polylib.cpp:17
Polylib::gammaF
double gammaF(const double)
Calculate the Gamma function , , for integer values and halves.
Definition: Polylib.cpp:1161
Polylib::Dgj
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
Definition: Polylib.cpp:546
Polylib::zwgrjm
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
Definition: Polylib.cpp:119
Polylib::JacZeros
void JacZeros(const int n, double *a, double *b, const double alpha, const double beta)
Zero and Weight determination through the eigenvalues and eigenvectors of a tridiagonal matrix from t...
Definition: Polylib.cpp:1336
Polylib::zwgrjp
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
Definition: Polylib.cpp:158
Polylib::laginterp
double laginterp(double z, int j, const double *zj, int np)
Definition: Polylib.cpp:42
Polylib::chri1
void chri1(int, double *, double *, double *, double *, double)
Given a weight function through the first n+1 coefficients a and b of its orthogonal polynomials thi...
Definition: Polylib.cpp:1608
Polylib::Dglj
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
Definition: Polylib.cpp:690
Polylib::Imgrjp
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at ...
Definition: Polylib.cpp:946
Polylib::zwglj
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
Definition: Polylib.cpp:195
Polylib::zwgj
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:91
Polylib::Imgj
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm.
Definition: Polylib.cpp:887
Polylib::RecCoeff
static void RecCoeff(const int, double *, double *, const double, const double)
The routine finds the recurrence coefficients a and b of the orthogonal polynomials.
Definition: Polylib.cpp:1367
Polylib::hgrjm
double hgrjm(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:787
Polylib::gammaFracGammaF
double gammaFracGammaF(const int, const double, const int, const double)
Calculate fraction of two Gamma functions, , for integer values and halves.
Definition: Polylib.cpp:1203
Polylib::TriQL
static void TriQL(const int, double *, double *, double **)
QL algorithm for symmetric tridiagonal matrix.
Definition: Polylib.cpp:1425
Polylib::Jacobz
static void Jacobz(const int n, double *z, const double alpha, const double beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e. .
Definition: Polylib.cpp:1281
Polylib::hgj
double hgj(const int i, const double z, const double *zgj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the ar...
Definition: Polylib.cpp:752
Polylib::Imgrjm
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at...
Definition: Polylib.cpp:916
Polylib::zwgk
void zwgk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:240
Polylib::Dgrjm
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
Definition: Polylib.cpp:589
Polylib::ImagBesselComp
std::complex< Nektar::NekDouble > ImagBesselComp(int n, std::complex< Nektar::NekDouble > y)
Calcualte the bessel function of the first kind with complex double input y. Taken from Numerical Rec...
Definition: Polylib.cpp:1651
Polylib::Dgrjp
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1.
Definition: Polylib.cpp:639
Polylib::hgrjp
double hgrjp(const int i, const double z, const double *zgrj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
Definition: Polylib.cpp:823
Polylib::zwrk
void zwrk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Radau-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:313
Polylib::hglj
double hglj(const int i, const double z, const double *zglj, const int np, const double alpha, const double beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zgrj ...
Definition: Polylib.cpp:858
Polylib::Imglj
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm.
Definition: Polylib.cpp:977
Polylib::jacobd
void jacobd(const int np, const double *z, double *polyd, const int n, const double alpha, const double beta)
Calculate the derivative of Jacobi polynomials.
Definition: Polylib.cpp:1134
Polylib::optdiff
double optdiff(double xl, double xr)
The following function is used to circumvent/reduce "Subtractive Cancellation" The expression 1/dz is...
Definition: Polylib.cpp:22
Polylib::jacobfd
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1034
Polylib::JKMatrix
void JKMatrix(int, double *, double *)
Calcualtes the Jacobi-kronrod matrix by determining the a and coefficients.
Definition: Polylib.cpp:1527
Polylib::zwlk
void zwlk(double *z, double *w, const int npt, const double alpha, const double beta)
Gauss-Lobatto-Kronrod-Jacobi zeros and weights.
Definition: Polylib.cpp:420
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:272
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:28 for Nektar++ by doxygen 1.9.1