Intrepid
http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/src/Shared/Intrepid_PolylibDef.hpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //
00005 //                           Intrepid Package
00006 //                 Copyright (2007) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00039 //                    Denis Ridzal  (dridzal@sandia.gov), or
00040 //                    Kara Peterson (kjpeter@sandia.gov)
00041 //
00042 // ************************************************************************
00043 // @HEADER
00044 */
00045 
00047 //
00048 // File: Intrepid_PolylibDef.hpp
00049 //
00050 // For more information, please see: http://www.nektar.info
00051 //
00052 // The MIT License
00053 //
00054 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
00055 // Department of Aeronautics, Imperial College London (UK), and Scientific
00056 // Computing and Imaging Institute, University of Utah (USA).
00057 //
00058 // License for the specific language governing rights and limitations under
00059 // Permission is hereby granted, free of charge, to any person obtaining a
00060 // copy of this software and associated documentation files (the "Software"),
00061 // to deal in the Software without restriction, including without limitation
00062 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
00063 // and/or sell copies of the Software, and to permit persons to whom the
00064 // Software is furnished to do so, subject to the following conditions:
00065 //
00066 // The above copyright notice and this permission notice shall be included
00067 // in all copies or substantial portions of the Software.
00068 //
00069 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
00070 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
00071 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
00072 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
00073 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
00074 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
00075 // DEALINGS IN THE SOFTWARE.
00076 //
00077 // Description:
00078 // This file is redistributed with the Intrepid package. It should be used
00079 // in accordance with the above MIT license, at the request of the authors.
00080 // This file is NOT covered by the usual Intrepid/Trilinos LGPL license.
00081 //
00082 // Origin: Nektar++ library, http://www.nektar.info, downloaded on
00083 //         March 10, 2009.
00084 //
00086 
00087 
00094 #ifdef _MSC_VER
00095 #include "winmath.h"
00096 #endif
00097 
00098 
00099 namespace Intrepid {
00100 
00102 #define INTREPID_POLYLIB_STOP 50 
00103 
00105 #define INTREPID_POLYLIB_POLYNOMIAL_DEFLATION 0
00106 
00107 #ifdef INTREPID_POLYLIB_POLYNOMIAL_DEFLATION
00108 
00109 #define jacobz(n,z,alpha,beta) Jacobz(n,z,alpha,beta)
00110 #else
00111 
00112 #define jacobz(n,z,alpha,beta) JacZeros(n,z,alpha,beta)
00113 #endif
00114 
00115 
00116 template <class Scalar>
00117 void IntrepidPolylib::zwgj (Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
00118   register int i;
00119   Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
00120 
00121   IntrepidPolylib::jacobz (np,z,alpha,beta);
00122   IntrepidPolylib::jacobd (np,z,w,np,alpha,beta);
00123 
00124   fac  = std::pow(two,apb + one)*IntrepidPolylib::gammaF(alpha + np + one)*IntrepidPolylib::gammaF(beta + np + one);
00125   fac /= IntrepidPolylib::gammaF((Scalar)(np + one))*IntrepidPolylib::gammaF(apb + np + one);
00126 
00127   for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));
00128 
00129   return;
00130 }
00131 
00132 
00133 template <class Scalar>
00134 void IntrepidPolylib::zwgrjm(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
00135 
00136   if(np == 1){
00137     z[0] = 0.0;
00138     w[0] = 2.0;
00139   }
00140   else{
00141     register int i;
00142     Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
00143 
00144     z[0] = -one;
00145     IntrepidPolylib::jacobz  (np-1,z+1,alpha,beta+1);
00146     IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
00147 
00148     fac  = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
00149     fac /= IntrepidPolylib::gammaF((Scalar)np)*(beta + np)*IntrepidPolylib::gammaF(apb + np + 1);
00150 
00151     for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
00152     w[0] *= (beta + one);
00153   }
00154 
00155   return;
00156 }
00157 
00158 
00159 template <class Scalar>
00160 void IntrepidPolylib::zwgrjp(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
00161 
00162   if(np == 1){
00163     z[0] = 0.0;
00164     w[0] = 2.0;
00165   }
00166   else{
00167     register int i;
00168     Scalar fac, one = 1.0, two = 2.0, apb = alpha + beta;
00169 
00170     IntrepidPolylib::jacobz  (np-1,z,alpha+1,beta);
00171     z[np-1] = one;
00172     IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
00173 
00174     fac  = std::pow(two,apb)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
00175     fac /= IntrepidPolylib::gammaF((Scalar)np)*(alpha + np)*IntrepidPolylib::gammaF(apb + np + 1);
00176 
00177     for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
00178     w[np-1] *= (alpha + one);
00179   }
00180 
00181   return;
00182 }
00183 
00184 
00185 template <class Scalar>
00186 void IntrepidPolylib::zwglj(Scalar *z, Scalar *w, const int np, const Scalar alpha, const Scalar beta){
00187 
00188   if( np == 1 ){
00189     z[0] = 0.0;
00190     w[0] = 2.0;
00191   }
00192   else{
00193     register int i;
00194     Scalar   fac, one = 1.0, apb = alpha + beta, two = 2.0;
00195 
00196     z[0]    = -one;
00197     z[np-1] =  one;
00198     IntrepidPolylib::jacobz  (np-2,z + 1,alpha + one,beta + one);
00199     IntrepidPolylib::jacobfd (np,z,w,(Scalar*)0,np-1,alpha,beta);
00200 
00201     fac  = std::pow(two,apb + 1)*IntrepidPolylib::gammaF(alpha + np)*IntrepidPolylib::gammaF(beta + np);
00202     fac /= (np-1)*IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha + beta + np + one);
00203 
00204     for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
00205     w[0]    *= (beta  + one);
00206     w[np-1] *= (alpha + one);
00207   }
00208 
00209   return;
00210 }
00211 
00212 
00213 template <class Scalar>
00214 void IntrepidPolylib::Dgj(Scalar *D,  const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
00215 {
00216 
00217     Scalar one = 1.0, two = 2.0;
00218 
00219     if (np <= 0){
00220         D[0] = 0.0;
00221     }
00222     else{
00223         register int i,j; 
00224         Scalar *pd;
00225 
00226         pd = (Scalar *)malloc(np*sizeof(Scalar));
00227         IntrepidPolylib::jacobd(np,z,pd,np,alpha,beta);
00228 
00229         for (i = 0; i < np; i++){
00230             for (j = 0; j < np; j++){
00231 
00232                 if (i != j) 
00233                     //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.
00234                     D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
00235                 else    
00236                     D[i*np+j] = (alpha - beta + (alpha + beta + two)*z[j])/
00237                     (two*(one - z[j]*z[j]));
00238             }
00239         }
00240         free(pd);
00241     }
00242     return;
00243 }
00244 
00245 
00246 template <class Scalar>
00247 void IntrepidPolylib::Dgrjm(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
00248 {
00249 
00250     if (np <= 0){
00251         D[0] = 0.0;
00252     }
00253     else{
00254         register int i, j; 
00255         Scalar   one = 1.0, two = 2.0;
00256         Scalar   *pd;
00257 
00258         pd  = (Scalar *)malloc(np*sizeof(Scalar));
00259 
00260         pd[0] = std::pow(-one,np-1)*IntrepidPolylib::gammaF((Scalar)np+beta+one);
00261         pd[0] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(beta+two);
00262         IntrepidPolylib::jacobd(np-1,z+1,pd+1,np-1,alpha,beta+1);
00263         for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);
00264 
00265         for (i = 0; i < np; i++) {
00266             for (j = 0; j < np; j++){
00267                 if (i != j) 
00268                     //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.
00269                     D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
00270                 else { 
00271                     if(j == 0)
00272                         D[i*np+j] = -(np + alpha + beta + one)*(np - one)/
00273                         (two*(beta + two));
00274                     else
00275                         D[i*np+j] = (alpha - beta + one + (alpha + beta + one)*z[j])/
00276                         (two*(one - z[j]*z[j]));
00277                 }
00278             }
00279         }
00280         free(pd);
00281     }
00282     return;
00283 }
00284 
00285 
00286 template <class Scalar>
00287 void IntrepidPolylib::Dgrjp(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
00288 {
00289 
00290     if (np <= 0){
00291         D[0] = 0.0;
00292     }
00293     else{
00294         register int i, j; 
00295         Scalar   one = 1.0, two = 2.0;
00296         Scalar   *pd;
00297 
00298         pd  = (Scalar *)malloc(np*sizeof(Scalar));
00299 
00300 
00301         IntrepidPolylib::jacobd(np-1,z,pd,np-1,alpha+1,beta);
00302         for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
00303         pd[np-1] = -IntrepidPolylib::gammaF((Scalar)np+alpha+one);
00304         pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np)*IntrepidPolylib::gammaF(alpha+two);
00305 
00306         for (i = 0; i < np; i++) {
00307             for (j = 0; j < np; j++){
00308                 if (i != j)
00309                     //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.
00310                     D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
00311                 else { 
00312                     if(j == np-1)
00313                         D[i*np+j] = (np + alpha + beta + one)*(np - one)/
00314                         (two*(alpha + two));
00315                     else
00316                         D[i*np+j] = (alpha - beta - one + (alpha + beta + one)*z[j])/
00317                         (two*(one - z[j]*z[j]));
00318                 }
00319             }
00320         }
00321         free(pd);
00322     }
00323     return;
00324 }
00325 
00326 
00327 template <class Scalar>
00328 void IntrepidPolylib::Dglj(Scalar *D, const Scalar *z, const int np, const Scalar alpha, const Scalar beta)
00329 {
00330 
00331     if (np <= 0){
00332         D[0] = 0.0;
00333     }
00334     else{
00335         register int i, j; 
00336         Scalar   one = 1.0, two = 2.0;
00337         Scalar   *pd;
00338 
00339         pd  = (Scalar *)malloc(np*sizeof(Scalar));
00340 
00341         pd[0]  = two*std::pow(-one,np)*IntrepidPolylib::gammaF((Scalar)np + beta);
00342         pd[0] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(beta + two);
00343         IntrepidPolylib::jacobd(np-2,z+1,pd+1,np-2,alpha+1,beta+1);
00344         for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
00345         pd[np-1]  = -two*IntrepidPolylib::gammaF((Scalar)np + alpha);
00346         pd[np-1] /= IntrepidPolylib::gammaF((Scalar)np - one)*IntrepidPolylib::gammaF(alpha + two);
00347 
00348         for (i = 0; i < np; i++) {
00349             for (j = 0; j < np; j++){
00350                 if (i != j)
00351                     //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.
00352                     D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
00353                 else { 
00354                     if (j == 0)
00355                         D[i*np+j] = (alpha - (np-1)*(np + alpha + beta))/(two*(beta+ two));
00356                     else if (j == np-1)
00357                         D[i*np+j] =-(beta - (np-1)*(np + alpha + beta))/(two*(alpha+ two));
00358                     else
00359                         D[i*np+j] = (alpha - beta + (alpha + beta)*z[j])/
00360                         (two*(one - z[j]*z[j]));
00361                 }
00362             }
00363         }
00364         free(pd);
00365     }
00366     return;
00367 }
00368 
00369 
00370 template <class Scalar>
00371 Scalar IntrepidPolylib::hgj (const int i, const Scalar z, const Scalar *zgj,
00372                              const int np, const Scalar alpha, const Scalar beta)
00373 {
00374 
00375     Scalar zi, dz, p, pd, h;
00376 
00377     zi  = *(zgj+i);
00378     dz  = z - zi;
00379     if (std::abs(dz) < INTREPID_TOL) return 1.0;
00380 
00381     IntrepidPolylib::jacobd (1, &zi, &pd , np, alpha, beta);
00382     IntrepidPolylib::jacobfd(1, &z , &p, (Scalar*)0 , np, alpha, beta);
00383     h = p/(pd*dz);
00384 
00385     return h;
00386 }
00387 
00388 
00389 template <class Scalar>
00390 Scalar IntrepidPolylib::hgrjm (const int i, const Scalar z, const Scalar *zgrj,
00391                                const int np, const Scalar alpha, const Scalar beta)
00392 {
00393 
00394     Scalar zi, dz, p, pd, h;
00395 
00396     zi  = *(zgrj+i);
00397     dz  = z - zi;
00398     if (std::abs(dz) < INTREPID_TOL) return 1.0;
00399 
00400     IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha, beta + 1);
00401     // need to use this routine in case zi = -1 or 1
00402     IntrepidPolylib::jacobd  (1, &zi, &pd, np-1, alpha, beta + 1);
00403     h = (1.0 + zi)*pd + p;
00404     IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0,  np-1, alpha, beta + 1);
00405     h = (1.0 + z )*p/(h*dz);
00406 
00407     return h;
00408 }
00409 
00410 
00411 template <class Scalar>
00412 Scalar IntrepidPolylib::hgrjp (const int i, const Scalar z, const Scalar *zgrj,
00413                                const int np, const Scalar alpha, const Scalar beta)
00414 {
00415 
00416     Scalar zi, dz, p, pd, h;
00417 
00418     zi  = *(zgrj+i);
00419     dz  = z - zi;
00420     if (std::abs(dz) < INTREPID_TOL) return 1.0;
00421 
00422     IntrepidPolylib::jacobfd (1, &zi, &p , (Scalar*)0, np-1, alpha+1, beta );
00423     // need to use this routine in case z = -1 or 1
00424     IntrepidPolylib::jacobd  (1, &zi, &pd, np-1, alpha+1, beta );
00425     h = (1.0 - zi)*pd - p;
00426     IntrepidPolylib::jacobfd (1, &z, &p, (Scalar*)0,  np-1, alpha+1, beta);
00427     h = (1.0 - z )*p/(h*dz);
00428 
00429     return h;
00430 }
00431 
00432 
00433 template <class Scalar>
00434 Scalar IntrepidPolylib::hglj (const int i, const Scalar z, const Scalar *zglj,
00435                               const int np, const Scalar alpha, const Scalar beta)
00436 {
00437     Scalar one = 1., two = 2.;
00438     Scalar zi, dz, p, pd, h;
00439 
00440     zi  = *(zglj+i);
00441     dz  = z - zi;
00442     if (std::abs(dz) < INTREPID_TOL) return 1.0;
00443 
00444     IntrepidPolylib::jacobfd(1, &zi, &p , (Scalar*)0, np-2, alpha + one, beta + one);
00445     // need to use this routine in case z = -1 or 1
00446     IntrepidPolylib::jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
00447     h = (one - zi*zi)*pd - two*zi*p;
00448     IntrepidPolylib::jacobfd(1, &z, &p, (Scalar*)0, np-2, alpha + one, beta + one);
00449     h = (one - z*z)*p/(h*dz);
00450 
00451     return h;
00452 }
00453 
00454 
00455 template <class Scalar>
00456 void IntrepidPolylib::Imgj(Scalar *im, const Scalar *zgj, const Scalar *zm, const int nz, 
00457                            const int mz, const Scalar alpha, const Scalar beta){
00458         Scalar zp;
00459         register int i, j;
00460 
00461         /* old Polylib code */
00462         for (i = 0; i < mz; ++i) {
00463             zp = zm[i];
00464             for (j = 0; j < nz; ++j) {
00465                 im[i*nz+j] = IntrepidPolylib::hgj(j, zp, zgj, nz, alpha, beta);
00466             }
00467         }
00468 
00469         /* original Nektar++ code: is this correct???
00470         for (i = 0; i < nz; ++i) {
00471             for (j = 0; j < mz; ++j)
00472             {
00473                 zp = zm[j];
00474                 im [i*mz+j] = IntrepidPolylib::hgj(i, zp, zgj, nz, alpha, beta);
00475             }
00476         }
00477         */
00478 
00479         return;
00480 }
00481 
00482 
00483 template <class Scalar>
00484 void IntrepidPolylib::Imgrjm(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz,
00485                              const int mz, const Scalar alpha, const Scalar beta)
00486 {
00487     Scalar zp;
00488     register int i, j;
00489 
00490     for (i = 0; i < mz; i++) {
00491       zp = zm[i];
00492       for (j = 0; j < nz; j++) {
00493         im[i*nz+j] = IntrepidPolylib::hgrjm(j, zp, zgrj, nz, alpha, beta);
00494       }
00495     }
00496 
00497     /* original Nektar++ code: is this correct???
00498     for (i = 0; i < nz; i++) {
00499         for (j = 0; j < mz; j++)
00500         {
00501             zp = zm[j];
00502             im [i*mz+j] = IntrepidPolylib::hgrjm(i, zp, zgrj, nz, alpha, beta);
00503         }
00504     }
00505     */
00506 
00507     return;
00508 }
00509 
00510 
00511 template <class Scalar>
00512 void IntrepidPolylib::Imgrjp(Scalar *im, const Scalar *zgrj, const Scalar *zm, const int nz, 
00513                              const int mz, const Scalar alpha, const Scalar beta)
00514 {
00515         Scalar zp;
00516         register int i, j;
00517 
00518         for (i = 0; i < mz; i++) {
00519           zp = zm[i];
00520           for (j = 0; j < nz; j++) {
00521             im [i*nz+j] = IntrepidPolylib::hgrjp(j, zp, zgrj, nz, alpha, beta);
00522           }
00523         }
00524 
00525         /* original Nektar++ code: is this correct?
00526         for (i = 0; i < nz; i++) {
00527             for (j = 0; j < mz; j++)
00528             {
00529                 zp = zm[j];
00530                 im [i*mz+j] = IntrepidPolylib::hgrjp(i, zp, zgrj, nz, alpha, beta);
00531             }
00532         }
00533         */
00534 
00535         return;
00536 }
00537 
00538 
00539 template <class Scalar>
00540 void IntrepidPolylib::Imglj(Scalar *im, const Scalar *zglj, const Scalar *zm, const int nz, 
00541                             const int mz, const Scalar alpha, const Scalar beta)
00542 {
00543     Scalar zp;
00544     register int i, j;
00545 
00546     for (i = 0; i < mz; i++) {
00547       zp = zm[i];
00548       for (j = 0; j < nz; j++) {
00549         im[i*nz+j] = IntrepidPolylib::hglj(j, zp, zglj, nz, alpha, beta);
00550       }
00551     }
00552 
00553     /* original Nektar++ code: is this correct?
00554     for (i = 0; i < nz; i++) {
00555         for (j = 0; j < mz; j++)
00556         {
00557             zp = zm[j];
00558             im[i*mz+j] = IntrepidPolylib::hglj(i, zp, zglj, nz, alpha, beta);
00559         }
00560     }
00561     */
00562 
00563     return;
00564 }
00565 
00566 
00567 template <class Scalar>
00568 void IntrepidPolylib::jacobfd(const int np, const Scalar *z, Scalar *poly_in, Scalar *polyd,
00569                               const int n, const Scalar alpha, const Scalar beta){
00570   register int i;
00571   Scalar  zero = 0.0, one = 1.0, two = 2.0;
00572 
00573   if(!np)
00574     return;
00575 
00576   if(n == 0){
00577     if(poly_in)
00578       for(i = 0; i < np; ++i)
00579         poly_in[i] = one;
00580     if(polyd)
00581       for(i = 0; i < np; ++i)
00582         polyd[i] = zero;
00583   }
00584   else if (n == 1){
00585     if(poly_in)
00586       for(i = 0; i < np; ++i)
00587         poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
00588     if(polyd)
00589       for(i = 0; i < np; ++i)
00590         polyd[i] = 0.5*(alpha + beta + two);
00591   }
00592   else{
00593     register int k;
00594     Scalar   a1,a2,a3,a4;
00595     Scalar   two = 2.0, apb = alpha + beta;
00596     Scalar   *poly, *polyn1,*polyn2;
00597 
00598     if(poly_in){ // switch for case of no poynomial function return
00599       polyn1 = (Scalar *)malloc(2*np*sizeof(Scalar));
00600       polyn2 = polyn1+np;
00601       poly   = poly_in;
00602     }
00603     else{
00604       polyn1 = (Scalar *)malloc(3*np*sizeof(Scalar));
00605       polyn2 = polyn1+np;
00606       poly   = polyn2+np;
00607     }
00608 
00609     for(i = 0; i < np; ++i){
00610       polyn2[i] = one;
00611       polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
00612     }
00613 
00614     for(k = 2; k <= n; ++k){
00615       a1 =  two*k*(k + apb)*(two*k + apb - two);
00616       a2 = (two*k + apb - one)*(alpha*alpha - beta*beta);
00617       a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb);
00618       a4 =  two*(k + alpha - one)*(k + beta - one)*(two*k + apb);
00619 
00620       a2 /= a1;
00621       a3 /= a1;
00622       a4 /= a1;
00623 
00624       for(i = 0; i < np; ++i){
00625         poly  [i] = (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i];
00626         polyn2[i] = polyn1[i];
00627         polyn1[i] = poly  [i];
00628       }
00629     }
00630 
00631     if(polyd){
00632       a1 = n*(alpha - beta);
00633       a2 = n*(two*n + alpha + beta);
00634       a3 = two*(n + alpha)*(n + beta);
00635       a4 = (two*n + alpha + beta);
00636       a1 /= a4;  a2 /= a4;   a3 /= a4;
00637 
00638       // note polyn2 points to polyn1 at end of poly iterations
00639       for(i = 0; i < np; ++i){
00640         polyd[i]  = (a1- a2*z[i])*poly[i] + a3*polyn2[i];
00641         polyd[i] /= (one - z[i]*z[i]);
00642       }
00643     }
00644 
00645     free(polyn1);
00646   }
00647 
00648   return;
00649 }
00650 
00651 
00652 template <class Scalar>
00653 void IntrepidPolylib::jacobd(const int np, const Scalar *z, Scalar *polyd, const int n,
00654                              const Scalar alpha, const Scalar beta)
00655 {
00656   register int i;
00657   Scalar one = 1.0;
00658   if(n == 0)
00659     for(i = 0; i < np; ++i) polyd[i] = 0.0;
00660   else{
00661     //jacobf(np,z,polyd,n-1,alpha+one,beta+one);
00662     IntrepidPolylib::jacobfd(np,z,polyd,(Scalar*)0,n-1,alpha+one,beta+one);
00663     for(i = 0; i < np; ++i) polyd[i] *= 0.5*(alpha + beta + (Scalar)n + one);
00664   }
00665   return;
00666 }
00667 
00668 
00669 template <class Scalar>
00670 void IntrepidPolylib::Jacobz(const int n, Scalar *z, const Scalar alpha, const Scalar beta){
00671     register int i,j,k;
00672     Scalar   dth = M_PI/(2.0*(Scalar)n);
00673     Scalar   poly,pder,rlast=0.0;
00674     Scalar   sum,delr,r;
00675     Scalar one = 1.0, two = 2.0;
00676 
00677     if(!n)
00678         return;
00679 
00680     for(k = 0; k < n; ++k){
00681         r = -std::cos((two*(Scalar)k + one) * dth);
00682         if(k) r = 0.5*(r + rlast);
00683 
00684         for(j = 1; j < INTREPID_POLYLIB_STOP; ++j){
00685             IntrepidPolylib::jacobfd(1,&r,&poly, &pder, n, alpha, beta);
00686 
00687             for(i = 0, sum = 0.0; i < k; ++i) sum += one/(r - z[i]);
00688 
00689             delr = -poly / (pder - sum * poly);
00690             r   += delr;
00691             if( std::abs(delr) < INTREPID_TOL ) break;
00692         }
00693         z[k]  = r;
00694         rlast = r;
00695     }
00696     return;
00697 }
00698 
00699 
00700 template <class Scalar>
00701 void IntrepidPolylib::JacZeros(const int n, Scalar *a, const Scalar alpha, const Scalar beta){
00702   int i;
00703   Scalar apb, apbi,a2b2;
00704   Scalar *b;
00705 
00706   if(!n)
00707     return;
00708 
00709   b = (Scalar *) malloc(n*sizeof(Scalar));
00710 
00711   // generate normalised terms
00712   apb  = alpha + beta;
00713   apbi = 2.0 + apb;
00714 
00715   b[n-1] = std::pow(2.0,apb+1.0)*IntrepidPolylib::gammaF(alpha+1.0)*IntrepidPolylib::gammaF(beta+1.0)/gammaF(apbi);
00716   a[0]   = (beta-alpha)/apbi;
00717   b[0]   = std::sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));
00718 
00719   a2b2 = beta*beta-alpha*alpha;
00720   for(i = 1; i < n-1; ++i){
00721     apbi = 2.0*(i+1) + apb;
00722     a[i] = a2b2/((apbi-2.0)*apbi);
00723     b[i] = std::sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
00724                      ((apbi*apbi-1)*apbi*apbi));
00725   }
00726 
00727   apbi   = 2.0*n + apb;
00728   //a[n-1] = a2b2/((apbi-2.0)*apbi); // THIS IS A BUG!!!
00729   if (n>1) a[n-1] = a2b2/((apbi-2.0)*apbi);
00730 
00731   // find eigenvalues
00732   IntrepidPolylib::TriQL(n, a, b);
00733 
00734   free(b);
00735   return;
00736 }
00737 
00738 
00739 template <class Scalar>
00740 void IntrepidPolylib::TriQL(const int n, Scalar *d,Scalar *e) {
00741   int m,l,iter,i,k;
00742   Scalar s,r,p,g,f,dd,c,b;
00743 
00744   for (l=0;l<n;l++) {
00745     iter=0;
00746     do {
00747       for (m=l;m<n-1;m++) {
00748         dd=std::abs(d[m])+std::abs(d[m+1]);
00749         if (std::abs(e[m])+dd == dd) break;
00750       }
00751       if (m != l) {
00752         if (iter++ == 50){
00753           TEUCHOS_TEST_FOR_EXCEPTION((1),
00754                              std::runtime_error,
00755                              ">>> ERROR (IntrepidPolylib): Too many iterations in TQLI.");
00756         }
00757         g=(d[l+1]-d[l])/(2.0*e[l]);
00758         r=std::sqrt((g*g)+1.0);
00759         //g=d[m]-d[l]+e[l]/(g+sign(r,g));
00760         g=d[m]-d[l]+e[l]/(g+((g)<0 ? -std::abs(r) : std::abs(r)));
00761         s=c=1.0;
00762         p=0.0;
00763         for (i=m-1;i>=l;i--) {
00764           f=s*e[i];
00765           b=c*e[i];
00766           if (std::abs(f) >= std::abs(g)) {
00767             c=g/f;
00768             r=std::sqrt((c*c)+1.0);
00769             e[i+1]=f*r;
00770             c *= (s=1.0/r);
00771           } else {
00772             s=f/g;
00773             r=std::sqrt((s*s)+1.0);
00774             e[i+1]=g*r;
00775             s *= (c=1.0/r);
00776           }
00777           g=d[i+1]-p;
00778           r=(d[i]-g)*s+2.0*c*b;
00779           p=s*r;
00780           d[i+1]=g+p;
00781           g=c*r-b;
00782         }
00783         d[l]=d[l]-p;
00784         e[l]=g;
00785         e[m]=0.0;
00786       }
00787     } while (m != l);
00788   }
00789 
00790   // order eigenvalues
00791   for(i = 0; i < n-1; ++i){
00792     k = i;
00793     p = d[i];
00794     for(l = i+1; l < n; ++l)
00795       if (d[l] < p) {
00796         k = l;
00797         p = d[l];
00798       }
00799     d[k] = d[i];
00800     d[i] = p;
00801   }
00802 }
00803 
00804 
00805 template <class Scalar>
00806 Scalar IntrepidPolylib::gammaF(const Scalar x){
00807   Scalar gamma = 1.0;
00808 
00809   if     (x == -0.5) gamma = -2.0*std::sqrt(M_PI);
00810   else if (!x) return gamma;
00811   else if ((x-(int)x) == 0.5){
00812     int n = (int) x;
00813     Scalar tmp = x;
00814 
00815     gamma = std::sqrt(M_PI);
00816     while(n--){
00817       tmp   -= 1.0;
00818       gamma *= tmp;
00819     }
00820   }
00821   else if ((x-(int)x) == 0.0){
00822     int n = (int) x;
00823     Scalar tmp = x;
00824 
00825     while(--n){
00826       tmp   -= 1.0;
00827       gamma *= tmp;
00828     }
00829   }
00830   else
00831     TEUCHOS_TEST_FOR_EXCEPTION((1),
00832                        std::invalid_argument,
00833                        ">>> ERROR (IntrepidPolylib): Argument is not of integer or half order.");
00834   return gamma;
00835 }
00836 
00837 } // end of namespace Intrepid
00838