Intrepid
Intrepid_PolylibDef.hpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ************************************************************************
4 //
5 // Intrepid Package
6 // Copyright (2007) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
39 // Denis Ridzal (dridzal@sandia.gov), or
40 // Kara Peterson (kjpeter@sandia.gov)
41 //
42 // ************************************************************************
43 // @HEADER
44 */
45 
47 //
48 // File: Intrepid_PolylibDef.hpp
49 //
50 // For more information, please see: http://www.nektar.info
51 //
52 // The MIT License
53 //
54 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
55 // Department of Aeronautics, Imperial College London (UK), and Scientific
56 // Computing and Imaging Institute, University of Utah (USA).
57 //
58 // License for the specific language governing rights and limitations under
59 // Permission is hereby granted, free of charge, to any person obtaining a
60 // copy of this software and associated documentation files (the "Software"),
61 // to deal in the Software without restriction, including without limitation
62 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
63 // and/or sell copies of the Software, and to permit persons to whom the
64 // Software is furnished to do so, subject to the following conditions:
65 //
66 // The above copyright notice and this permission notice shall be included
67 // in all copies or substantial portions of the Software.
68 //
69 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
70 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
71 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
72 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
73 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
74 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
75 // DEALINGS IN THE SOFTWARE.
76 //
77 // Description:
78 // This file is redistributed with the Intrepid package. It should be used
79 // in accordance with the above MIT license, at the request of the authors.
80 // This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
81 //
82 // Origin: Nektar++ library, http://www.nektar.info, downloaded on
83 // March 10, 2009.
84 //
86 
87 
95 namespace Intrepid {
96 
98 #define INTREPID_POLYLIB_STOP 50
99 
101 #define INTREPID_POLYLIB_POLYNOMIAL_DEFLATION 0
102 
103 #ifdef INTREPID_POLYLIB_POLYNOMIAL_DEFLATION
104 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
106 #else
107 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
109 #endif
110 
111 
112 template <class Scalar>
113 void IntrepidPolylib::zwgj (Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
114  int i;
115  Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
116 
117  IntrepidPolylib::jacobz (np,z,alpha,beta);
118  IntrepidPolylib::jacobd (np,z,w,np,alpha,beta);
119 
120  fac = std::pow(two,apb + one)*IntrepidPolylib::gammaF(alpha + np + one)*IntrepidPolylib::gammaF(beta + np + one);
121  fac /= IntrepidPolylib::gammaF((Scalar)(np + one))*IntrepidPolylib::gammaF(apb + np + one);
122 
123  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
124 
125  return;
126 }
127 
128 
129 template <class Scalar>
130 void IntrepidPolylib::zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
131 
132  if(np == 1){
133  z[0] = 0.0;
134  w[0] = 2.0;
135  }
136  else{
137  int i;
138  Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
139 
140  z[0] = -one;
141  IntrepidPolylib::jacobz (np-1,z+1,alpha,beta+1);
142  IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
143 
144  fac = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
145  fac /= IntrepidPolylib::gammaF((Scalar)np)*(beta + np)*IntrepidPolylib::gammaF(apb + np + 1);
146 
147  for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
148  w[0] *= (beta + one);
149  }
150 
151  return;
152 }
153 
154 
155 template <class Scalar>
156 void IntrepidPolylib::zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
157 
158  if(np == 1){
159  z[0] = 0.0;
160  w[0] = 2.0;
161  }
162  else{
163  int i;
164  Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
165 
166  IntrepidPolylib::jacobz (np-1,z,alpha+1,beta);
167  z[np-1] = one;
168  IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
169 
170  fac = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
171  fac /= IntrepidPolylib::gammaF((Scalar)np)*(alpha + np)*IntrepidPolylib::gammaF(apb + np + 1);
172 
173  for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
174  w[np-1] *= (alpha + one);
175  }
176 
177  return;
178 }
179 
180 
181 template <class Scalar>
182 void IntrepidPolylib::zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
183 
184  if( np == 1 ){
185  z[0] = 0.0;
186  w[0] = 2.0;
187  }
188  else{
189  int i;
190  Scalar fac, one = 1.0, apb = alpha + beta, two = 2.0;
191 
192  z[0] = -one;
193  z[np-1] = one;
194  IntrepidPolylib::jacobz (np-2,z + 1,alpha + one,beta + one);
195  IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
196 
197  fac = std::pow(two,apb + 1)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
198  fac /= (np-1)*IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha + beta + np + one);
199 
200  for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
201  w[0] *= (beta + one);
202  w[np-1] *= (alpha + one);
203  }
204 
205  return;
206 }
207 
208 
209 template <class Scalar>
210 void IntrepidPolylib::Dgj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
211 {
212 
213  Scalar one = 1.0, two = 2.0;
214 
215  if (np <= 0){
216  D[0] = 0.0;
217  }
218  else{
219  int i,j;
220  Scalar *pd;
221 
222  pd = (Scalar *)malloc(np*sizeof(Scalar));
223  IntrepidPolylib::jacobd(np,z,pd,np,alpha,beta);
224 
225  for (i = 0; i < np; i++){
226  for (j = 0; j < np; j++){
227 
228  if (i != j)
229  //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
230  D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
231  else
232  D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
233  (two*(one - z[j]*z[j]));
234  }
235  }
236  free(pd);
237  }
238  return;
239 }
240 
241 
242 template <class Scalar>
243 void IntrepidPolylib::Dgrjm(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
244 {
245 
246  if (np <= 0){
247  D[0] = 0.0;
248  }
249  else{
250  int i, j;
251  Scalar one = 1.0, two = 2.0;
252  Scalar *pd;
253 
254  pd = (Scalar *)malloc(np*sizeof(Scalar));
255 
256  pd[0] = std::pow(-one,np-1)*IntrepidPolylib::gammaF((Scalar)np+beta+one);
257  pd[0] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(beta+two);
258  IntrepidPolylib::jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
259  for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
260 
261  for (i = 0; i < np; i++) {
262  for (j = 0; j < np; j++){
263  if (i != j)
264  //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
265  D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
266  else {
267  if(j == 0)
268  D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
269  (two*(beta + two));
270  else
271  D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
272  (two*(one - z[j]*z[j]));
273  }
274  }
275  }
276  free(pd);
277  }
278  return;
279 }
280 
281 
282 template <class Scalar>
283 void IntrepidPolylib::Dgrjp(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
284 {
285 
286  if (np <= 0){
287  D[0] = 0.0;
288  }
289  else{
290  int i, j;
291  Scalar one = 1.0, two = 2.0;
292  Scalar *pd;
293 
294  pd = (Scalar *)malloc(np*sizeof(Scalar));
295 
296 
297  IntrepidPolylib::jacobd(np-1,z,pd,np-1,alpha+1,beta);
298  for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
299  pd[np-1] = -IntrepidPolylib::gammaF((Scalar)np+alpha+one);
300  pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha+two);
301 
302  for (i = 0; i < np; i++) {
303  for (j = 0; j < np; j++){
304  if (i != j)
305  //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
306  D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
307  else {
308  if(j == np-1)
309  D[i*np+j] = (np + alpha + beta + one)*(np - one)/
310  (two*(alpha + two));
311  else
312  D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
313  (two*(one - z[j]*z[j]));
314  }
315  }
316  }
317  free(pd);
318  }
319  return;
320 }
321 
322 
323 template <class Scalar>
324 void IntrepidPolylib::Dglj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
325 {
326 
327  if (np <= 0){
328  D[0] = 0.0;
329  }
330  else{
331  int i, j;
332  Scalar one = 1.0, two = 2.0;
333  Scalar *pd;
334 
335  pd = (Scalar *)malloc(np*sizeof(Scalar));
336 
337  pd[0] = two*std::pow(-one,np)*IntrepidPolylib::gammaF((Scalar)np + beta);
338  pd[0] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(beta + two);
339  IntrepidPolylib::jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
340  for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
341  pd[np-1] = -two*IntrepidPolylib::gammaF((Scalar)np + alpha);
342  pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(alpha + two);
343 
344  for (i = 0; i < np; i++) {
345  for (j = 0; j < np; j++){
346  if (i != j)
347  //D[i*np+j] = pd[j]/(pd[i]*(z[j]-z[i])); <--- This is either a bug, or the derivative matrix is not defined consistently.
348  D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
349  else {
350  if (j == 0)
351  D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
352  else if (j == np-1)
353  D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
354  else
355  D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
356  (two*(one - z[j]*z[j]));
357  }
358  }
359  }
360  free(pd);
361  }
362  return;
363 }
364 
365 
366 template <class Scalar>
367 Scalar IntrepidPolylib::hgj (const int i, const Scalar z, const Scalar *zgj,
368  const int np, const Scalar alpha, const Scalar beta)
369 {
370 
371  Scalar zi, dz, p, pd, h;
372 
373  zi = *(zgj+i);
374  dz = z - zi;
375  if (std::abs(dz) < INTREPID_TOL) return 1.0;
376 
377  IntrepidPolylib::jacobd (1, &zi, &pd , np, alpha, beta);
378  IntrepidPolylib::jacobfd(1, &z , &p, (Scalar*)0 , np, alpha, beta);
379  h = p/(pd*dz);
380 
381  return h;
382 }
383 
384 
385 template <class Scalar>
386 Scalar IntrepidPolylib::hgrjm (const int i, const Scalar z, const Scalar *zgrj,
387  const int np, const Scalar alpha, const Scalar beta)
388 {
389 
390  Scalar zi, dz, p, pd, h;
391 
392  zi = *(zgrj+i);
393  dz = z - zi;
394  if (std::abs(dz) < INTREPID_TOL) return 1.0;
395 
396  IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha, beta + 1);
397  // need to use this routine in case zi = -1 or 1
398  IntrepidPolylib::jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
399  h = (1.0 + zi)*pd + p;
400  IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0, np-1, alpha, beta + 1);
401  h = (1.0 + z )*p/(h*dz);
402 
403  return h;
404 }
405 
406 
407 template <class Scalar>
408 Scalar IntrepidPolylib::hgrjp (const int i, const Scalar z, const Scalar *zgrj,
409  const int np, const Scalar alpha, const Scalar beta)
410 {
411 
412  Scalar zi, dz, p, pd, h;
413 
414  zi = *(zgrj+i);
415  dz = z - zi;
416  if (std::abs(dz) < INTREPID_TOL) return 1.0;
417 
418  IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha+1, beta );
419  // need to use this routine in case z = -1 or 1
420  IntrepidPolylib::jacobd (1, &zi, &pd, np-1, alpha+1, beta );
421  h = (1.0 - zi)*pd - p;
422  IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0, np-1, alpha+1, beta);
423  h = (1.0 - z )*p/(h*dz);
424 
425  return h;
426 }
427 
428 
429 template <class Scalar>
430 Scalar IntrepidPolylib::hglj (const int i, const Scalar z, const Scalar *zglj,
431  const int np, const Scalar alpha, const Scalar beta)
432 {
433  Scalar one = 1., two = 2.;
434  Scalar zi, dz, p, pd, h;
435 
436  zi = *(zglj+i);
437  dz = z - zi;
438  if (std::abs(dz) < INTREPID_TOL) return 1.0;
439 
440  IntrepidPolylib::jacobfd(1, &zi, &p , (Scalar*)0, np-2, alpha + one, beta + one);
441  // need to use this routine in case z = -1 or 1
442  IntrepidPolylib::jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
443  h = (one - zi*zi)*pd - two*zi*p;
444  IntrepidPolylib::jacobfd(1, &z, &p, (Scalar*)0, np-2, alpha + one, beta + one);
445  h = (one - z*z)*p/(h*dz);
446 
447  return h;
448 }
449 
450 
451 template <class Scalar>
452 void IntrepidPolylib::Imgj(Scalar *im, const Scalar *zgj, const Scalar *zm, const int nz,
453  const int mz, const Scalar alpha, const Scalar beta){
454  Scalar zp;
455  int i, j;
456 
457  /* old Polylib code */
458  for (i = 0; i < mz; ++i) {
459  zp = zm[i];
460  for (j = 0; j < nz; ++j) {
461  im[i*nz+j] = IntrepidPolylib::hgj(j, zp, zgj, nz, alpha, beta);
462  }
463  }
464 
465  /* original Nektar++ code: is this correct???
466  for (i = 0; i < nz; ++i) {
467  for (j = 0; j < mz; ++j)
468  {
469  zp = zm[j];
470  im [i*mz+j] = IntrepidPolylib::hgj(i, zp, zgj, nz, alpha, beta);
471  }
472  }
473  */
474 
475  return;
476 }
477 
478 
479 template <class Scalar>
480 void IntrepidPolylib::Imgrjm(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz,
481  const int mz, const Scalar alpha, const Scalar beta)
482 {
483  Scalar zp;
484  int i, j;
485 
486  for (i = 0; i < mz; i++) {
487  zp = zm[i];
488  for (j = 0; j < nz; j++) {
489  im[i*nz+j] = IntrepidPolylib::hgrjm(j, zp, zgrj, nz, alpha, beta);
490  }
491  }
492 
493  /* original Nektar++ code: is this correct???
494  for (i = 0; i < nz; i++) {
495  for (j = 0; j < mz; j++)
496  {
497  zp = zm[j];
498  im [i*mz+j] = IntrepidPolylib::hgrjm(i, zp, zgrj, nz, alpha, beta);
499  }
500  }
501  */
502 
503  return;
504 }
505 
506 
507 template <class Scalar>
508 void IntrepidPolylib::Imgrjp(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz,
509  const int mz, const Scalar alpha, const Scalar beta)
510 {
511  Scalar zp;
512  int i, j;
513 
514  for (i = 0; i < mz; i++) {
515  zp = zm[i];
516  for (j = 0; j < nz; j++) {
517  im [i*nz+j] = IntrepidPolylib::hgrjp(j, zp, zgrj, nz, alpha, beta);
518  }
519  }
520 
521  /* original Nektar++ code: is this correct?
522  for (i = 0; i < nz; i++) {
523  for (j = 0; j < mz; j++)
524  {
525  zp = zm[j];
526  im [i*mz+j] = IntrepidPolylib::hgrjp(i, zp, zgrj, nz, alpha, beta);
527  }
528  }
529  */
530 
531  return;
532 }
533 
534 
535 template <class Scalar>
536 void IntrepidPolylib::Imglj(Scalar *im, const Scalar *zglj, const Scalar *zm, const int nz,
537  const int mz, const Scalar alpha, const Scalar beta)
538 {
539  Scalar zp;
540  int i, j;
541 
542  for (i = 0; i < mz; i++) {
543  zp = zm[i];
544  for (j = 0; j < nz; j++) {
545  im[i*nz+j] = IntrepidPolylib::hglj(j, zp, zglj, nz, alpha, beta);
546  }
547  }
548 
549  /* original Nektar++ code: is this correct?
550  for (i = 0; i < nz; i++) {
551  for (j = 0; j < mz; j++)
552  {
553  zp = zm[j];
554  im[i*mz+j] = IntrepidPolylib::hglj(i, zp, zglj, nz, alpha, beta);
555  }
556  }
557  */
558 
559  return;
560 }
561 
562 
563 template <class Scalar>
564 void
566 jacobfd (const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd,
567  const int n, const Scalar alpha, const Scalar beta)
568 {
569  const Scalar zero = 0.0, one = 1.0, two = 2.0;
570 
571  if (! np) {
572  return;
573  }
574 
575  if (n == 0) {
576  if (poly_in) {
577  for (int i = 0; i < np; ++i) {
578  poly_in[i] = one;
579  }
580  }
581  if (polyd) {
582  for (int i = 0; i < np; ++i) {
583  polyd[i] = zero;
584  }
585  }
586  }
587  else if (n == 1) {
588  if (poly_in) {
589  for (int i = 0; i < np; ++i) {
590  poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
591  }
592  }
593  if (polyd) {
594  for (int i = 0; i < np; ++i) {
595  polyd[i] = 0.5*(alpha + beta + two);
596  }
597  }
598  }
599  else {
600  Scalar a1,a2,a3,a4;
601  Scalar apb = alpha + beta;
602  Scalar *poly, *polyn1,*polyn2;
603 
604  if (poly_in) { // switch for case of no poynomial function return
605  polyn1 = (Scalar *)malloc(2*np*sizeof(Scalar));
606  polyn2 = polyn1+np;
607  poly = poly_in;
608  }
609  else{
610  polyn1 = (Scalar *)malloc(3*np*sizeof(Scalar));
611  polyn2 = polyn1+np;
612  poly = polyn2+np;
613  }
614 
615  for (int i = 0; i < np; ++i) {
616  polyn2[i] = one;
617  polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
618  }
619 
620  for (int k = 2; k <= n; ++k) {
621  a1 = two*k*(k + apb)*(two*k + apb - two);
622  a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
623  a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
624  a4 = two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
625 
626  a2 /= a1;
627  a3 /= a1;
628  a4 /= a1;
629 
630  for (int i = 0; i < np; ++i) {
631  poly [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
632  polyn2[i] = polyn1[i];
633  polyn1[i] = poly [i];
634  }
635  }
636 
637  if (polyd) {
638  a1 = n*(alpha - beta);
639  a2 = n*(two*n + alpha + beta);
640  a3 = two*(n + alpha)*(n + beta);
641  a4 = (two*n + alpha + beta);
642  a1 /= a4; a2 /= a4; a3 /= a4;
643 
644  // note polyn2 points to polyn1 at end of poly iterations
645  for (int i = 0; i < np; ++i) {
646  polyd[i] = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
647  polyd[i] /= (one - z[i]*z[i]);
648  }
649  }
650 
651  free(polyn1);
652  }
653 }
654 
655 
656 template <class Scalar>
657 void IntrepidPolylib::jacobd(const int np, const Scalar *z, Scalar *polyd, const int n,
658  const Scalar alpha, const Scalar beta)
659 {
660  int i;
661  Scalar one = 1.0;
662  if(n == 0)
663  for(i = 0; i < np; ++i) polyd[i] = 0.0;
664  else{
665  //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
666  IntrepidPolylib::jacobfd(np,z,polyd,(Scalar*)0,n-1,alpha+one,beta+one);
667  for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (Scalar)n + one);
668  }
669  return;
670 }
671 
672 
673 template <class Scalar>
674 void IntrepidPolylib::Jacobz(const int n, Scalar *z, const Scalar alpha, const Scalar beta){
675  int i,j,k;
676  Scalar dth = M_PI/(2.0*(Scalar)n);
677  Scalar poly,pder,rlast=0.0;
678  Scalar sum,delr,r;
679  Scalar one = 1.0, two = 2.0;
680 
681  if(!n)
682  return;
683 
684  for(k = 0; k < n; ++k){
685  r = -std::cos((two*(Scalar)k + one) * dth);
686  if(k) r = 0.5*(r + rlast);
687 
688  for(j = 1; j < INTREPID_POLYLIB_STOP; ++j){
689  IntrepidPolylib::jacobfd(1,&r,&poly, &pder, n, alpha, beta);
690 
691  for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
692 
693  delr = -poly / (pder - sum * poly);
694  r += delr;
695  if( std::abs(delr) < INTREPID_TOL ) break;
696  }
697  z[k] = r;
698  rlast = r;
699  }
700  return;
701 }
702 
703 
704 template <class Scalar>
705 void IntrepidPolylib::JacZeros(const int n, Scalar *a, const Scalar alpha, const Scalar beta){
706  int i;
707  Scalar apb, apbi,a2b2;
708  Scalar *b;
709 
710  if(!n)
711  return;
712 
713  b = (Scalar *) malloc(n*sizeof(Scalar));
714 
715  // generate normalised terms
716  apb = alpha + beta;
717  apbi = 2.0 + apb;
718 
719  b[n-1] = std::pow(2.0,apb+1.0)*IntrepidPolylib::gammaF(alpha+1.0)*IntrepidPolylib::gammaF(beta+1.0)/gammaF(apbi);
720  a[0] = (beta-alpha)/apbi;
721  b[0] = std::sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
722 
723  a2b2 = beta*beta-alpha*alpha;
724  for(i = 1; i < n-1; ++i){
725  apbi = 2.0*(i+1) + apb;
726  a[i] = a2b2/((apbi-2.0)*apbi);
727  b[i] = std::sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
728  ((apbi*apbi-1)*apbi*apbi));
729  }
730 
731  apbi = 2.0*n + apb;
732  //a[n-1] = a2b2/((apbi-2.0)*apbi); // THIS IS A BUG!!!
733  if (n>1) a[n-1] = a2b2/((apbi-2.0)*apbi);
734 
735  // find eigenvalues
736  IntrepidPolylib::TriQL(n, a, b);
737 
738  free(b);
739  return;
740 }
741 
742 
743 template <class Scalar>
744 void IntrepidPolylib::TriQL(const int n, Scalar *d,Scalar *e) {
745  int m,l,iter,i,k;
746  Scalar s,r,p,g,f,dd,c,b;
747 
748  for (l=0;l<n;l++) {
749  iter=0;
750  do {
751  for (m=l;m<n-1;m++) {
752  dd=std::abs(d[m])+std::abs(d[m+1]);
753  if (std::abs(e[m])+dd == dd) break;
754  }
755  if (m != l) {
756  if (iter++ == 50){
757  TEUCHOS_TEST_FOR_EXCEPTION((1),
758  std::runtime_error,
759  ">>> ERROR (IntrepidPolylib): Too many iterations in TQLI.");
760  }
761  g=(d[l+1]-d[l])/(2.0*e[l]);
762  r=std::sqrt((g*g)+1.0);
763  //g=d[m]-d[l]+e[l]/(g+sign(r,g));
764  g=d[m]-d[l]+e[l]/(g+((g)<0 ? -std::abs(r) : std::abs(r)));
765  s=c=1.0;
766  p=0.0;
767  for (i=m-1;i>=l;i--) {
768  f=s*e[i];
769  b=c*e[i];
770  if (std::abs(f) >= std::abs(g)) {
771  c=g/f;
772  r=std::sqrt((c*c)+1.0);
773  e[i+1]=f*r;
774  c *= (s=1.0/r);
775  } else {
776  s=f/g;
777  r=std::sqrt((s*s)+1.0);
778  e[i+1]=g*r;
779  s *= (c=1.0/r);
780  }
781  g=d[i+1]-p;
782  r=(d[i]-g)*s+2.0*c*b;
783  p=s*r;
784  d[i+1]=g+p;
785  g=c*r-b;
786  }
787  d[l]=d[l]-p;
788  e[l]=g;
789  e[m]=0.0;
790  }
791  } while (m != l);
792  }
793 
794  // order eigenvalues
795  for(i = 0; i < n-1; ++i){
796  k = i;
797  p = d[i];
798  for(l = i+1; l < n; ++l)
799  if (d[l] < p) {
800  k = l;
801  p = d[l];
802  }
803  d[k] = d[i];
804  d[i] = p;
805  }
806 }
807 
808 
809 template <class Scalar>
810 Scalar IntrepidPolylib::gammaF(const Scalar x){
811  Scalar gamma = 1.0;
812 
813  if (x == -0.5) gamma = -2.0*std::sqrt(M_PI);
814  else if (!x) return gamma;
815  else if ((x-(int)x) == 0.5){
816  int n = (int) x;
817  Scalar tmp = x;
818 
819  gamma = std::sqrt(M_PI);
820  while(n--){
821  tmp -= 1.0;
822  gamma *= tmp;
823  }
824  }
825  else if ((x-(int)x) == 0.0){
826  int n = (int) x;
827  Scalar tmp = x;
828 
829  while(--n){
830  tmp -= 1.0;
831  gamma *= tmp;
832  }
833  }
834  else
835  TEUCHOS_TEST_FOR_EXCEPTION((1),
836  std::invalid_argument,
837  ">>> ERROR (IntrepidPolylib): Argument is not of integer or half order.");
838  return gamma;
839 }
840 
841 } // end of namespace Intrepid
842 
843 
844 #if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
845 #ifdef __GNUC__
846 #warning "The Intrepid package is deprecated"
847 #endif
848 #endif
849 
static void Dgrjm(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Radau-Jacobi zeros with a z...
static void Imgj(Scalar *im, const Scalar *zgj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Jacobi points to an arbitrary distribution at points zm...
static Scalar hgrjp(const int i, const Scalar z, const Scalar *zgrj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
static void zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
static void TriQL(const int n, Scalar *d, Scalar *e)
QL algorithm for symmetric tridiagonal matrix.
#define INTREPID_POLYLIB_STOP
Maximum number of iterations in polynomial defalation routine Jacobz.
static Scalar hglj(const int i, const Scalar z, const Scalar *zglj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Lobatto-Jacobi points zglj ...
static void zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
static Scalar hgj(const int i, const Scalar z, const Scalar *zgj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Jacobi points zgj at the ar...
static void zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
static Scalar hgrjm(const int i, const Scalar z, const Scalar *zgrj, const int np, const Scalar alpha, const Scalar beta)
Compute the value of the i th Lagrangian interpolant through the np Gauss-Radau-Jacobi points zgrj at...
static void JacZeros(const int n, Scalar *a, const Scalar alpha, const Scalar beta)
Zero determination through the eigenvalues of a tridiagonal matrix from the three term recursion rela...
static void jacobd(const int np, const Scalar *z, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Calculate the derivative of Jacobi polynomials.
static void Dgrjp(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix associated with the Gauss-Radau-Jacobi zeros with a zero at z=1...
static void Imgrjm(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=-1) to an arbitrary distrubtion at...
static void jacobfd(const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd, const int n, const Scalar alpha, const Scalar beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
static void Jacobz(const int n, Scalar *z, const Scalar alpha, const Scalar beta)
Calculate the n zeros, z, of the Jacobi polynomial, i.e. .
static void Dgj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix and its transpose associated with the Gauss-Jacobi zeros.
static void Dglj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
Compute the Derivative Matrix associated with the Gauss-Lobatto-Jacobi zeros.
static void zwgj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta)
Gauss-Jacobi zeros and weights.
static void Imglj(Scalar *im, const Scalar *zglj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points to an arbitrary distrubtion at points zm...
static Scalar gammaF(const Scalar x)
Calculate the Gamma function , , for integer values x and halves.
#define jacobz(n, z, alpha, beta)
zero determination using Newton iteration with polynomial deflation
static void Imgrjp(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz, const int mz, const Scalar alpha, const Scalar beta)
Interpolation Operator from Gauss-Radau-Jacobi points (including z=1) to an arbitrary distrubtion at ...