00001 #include "SundanceFeketeQuadQuadrature.hpp"
00002 #include "SundanceOut.hpp"
00003 #include "SundancePoint.hpp"
00004 #include "SundanceGaussLobatto1D.hpp"
00005 #include "PlayaTabs.hpp"
00006
00007 using namespace Sundance;
00008 using namespace Teuchos;
00009
00010 extern "C"
00011 {
00012
00013
00014 void dgetrf_(const int* M, const int* N, double* A, const int* lda,
00015 const int* iPiv, int* info);
00016
00017
00018 void dgetri_(const int* n, double* a, const int* lda, const int* iPiv,
00019 double* work, const int* lwork, int* info);
00020 }
00021
00022 void FeketeQuadQuadrature::getPoints(int order, Array<double>& wgt, Array<
00023 double>& x, Array<double>& y)
00024 {
00025 int p = order + 3;
00026 p = p + (p % 2);
00027 int nNodes = p / 2;
00028 GaussLobatto1D rule(nNodes, 0.0, 1.0);
00029 Array<double> s = rule.nodes();
00030 Array<double> t = s;
00031 Array<double> w = rule.weights();
00032 int n = rule.nPoints();
00033
00034 wgt.resize(n * n);
00035 x.resize(n * n);
00036 y.resize(n * n);
00037
00038 int k = 0;
00039 for (int i = 0; i < n; i++)
00040 {
00041 double p = s[i];
00042 for (int j = 0; j < n; j++, k++)
00043 {
00044 double q = t[j];
00045 x[k] = p;
00046 y[k] = q;
00047 wgt[k] = w[i] * w[j];
00048 }
00049 }
00050 }
00051
00052 void FeketeQuadQuadrature::computeBasisCoeffs(const int order, Array<double>& basisCoeffs)
00053 {
00054
00055 Array<Point> feketePts;
00056
00057 Array<double> x;
00058 Array<double> y;
00059 Array<double> w;
00060 getPoints(order, w, x, y);
00061 int nFeketePts = w.length();
00062 feketePts.resize(nFeketePts);
00063 for (int i = 0; i < nFeketePts; i++)
00064 feketePts[i] = Point(x[i], y[i]);
00065
00066
00067 basisCoeffs.resize(nFeketePts * nFeketePts);
00068
00069
00070
00071 for (int n = 0; n < nFeketePts; n++)
00072 {
00073
00074
00075 double* start = &(basisCoeffs[n * nFeketePts]);
00076 evalPolynomials(nFeketePts, feketePts[n][0], feketePts[n][1], start);
00077 }
00078
00079
00080
00081 int lapack_err = 0;
00082 Array<int> pivot;
00083 pivot.resize(nFeketePts);
00084 Array<double> work;
00085 work.resize(1);
00086 int lwork = -1;
00087
00088
00089 ::dgetrf_(&nFeketePts, &nFeketePts, &(basisCoeffs[0]), &nFeketePts,
00090 &(pivot[0]), &lapack_err);
00091
00092 TEUCHOS_TEST_FOR_EXCEPTION(
00093 lapack_err != 0,
00094 std::runtime_error,
00095 "FeketeQuadQuadrature::computeBasisCoeffs(): factorization of generalized Vandermonde matrix failed");
00096
00097
00098 ::dgetri_(&nFeketePts, &(basisCoeffs[0]), &nFeketePts, &(pivot[0]),
00099 &(work[0]), &lwork, &lapack_err);
00100 lwork = (int) work[0];
00101 work.resize(lwork);
00102 ::dgetri_(&nFeketePts, &(basisCoeffs[0]), &nFeketePts, &(pivot[0]),
00103 &(work[0]), &lwork, &lapack_err);
00104
00105 TEUCHOS_TEST_FOR_EXCEPTION(
00106 lapack_err != 0,
00107 std::runtime_error,
00108 "FeketeQuadQuadrature::computeBasisCoeffs(): inversion of generalized Vandermonde matrix failed");
00109 }
00110
00111 void FeketeQuadQuadrature::evalPolynomials(int nPts, double x,
00112 double y, double* resultPtr)
00113 {
00114
00115 int i = 0;
00116 int order = (int) sqrt((double) nPts);
00117
00118 double xLegendre = 1.;
00119 double xLegendre1 = 0.;
00120
00121 for (int k = 0; k < order; k++)
00122 {
00123 double yLegendre = 1.;
00124 double yLegendre1 = 0.;
00125 double normxLegendre = sqrt( 2. * k + 1 ) * xLegendre;
00126
00127
00128 for (int l = 0; l < order ; l++)
00129 {
00130 double normyLegendre = sqrt( 2. * l + 1 ) * yLegendre;
00131
00132
00133 resultPtr[i++] = normxLegendre * normyLegendre;
00134
00135 double yLegendre2 = yLegendre1;
00136 yLegendre1 = yLegendre;
00137 yLegendre = ( (2 * l + 1) * (2 * y - 1) * yLegendre1 - l * yLegendre2 ) / ( l + 1 );
00138 }
00139
00140 double xLegendre2 = xLegendre1;
00141 xLegendre1 = xLegendre;
00142 xLegendre = ( (2 * k + 1) * (2 * x - 1) * xLegendre1 - k * xLegendre2 ) / ( k + 1 );
00143 }
00144 }
00145
00146 bool FeketeQuadQuadrature::test(int p)
00147 {
00148 Array<double> w;
00149 Array<double> x;
00150 Array<double> y;
00151
00152 getPoints(p, w, x, y);
00153 bool pass = true;
00154 for (int a = 0; a <= p; a++)
00155 {
00156 int bMax = p - a;
00157 for (int b = 0; b <= bMax; b++)
00158 {
00159 double sum = 0.0;
00160 for (int q = 0; q < w.length(); q++)
00161 {
00162 sum += w[q] * pow(x[q], (double) a) * pow(y[q], (double) b);
00163 }
00164 double err = fabs(sum - exact(a, b));
00165 bool localPass = err < 1.0e-14;
00166 pass = pass && localPass;
00167 if (!localPass)
00168 {
00169 fprintf(stderr,
00170 "order=%d m (%d, %d) q=%22.15g exact=%22.15g\n", p, a,
00171 b, sum, exact(a, b));
00172 std::cerr << "error = " << err << std::endl;
00173 }
00174 }
00175 }
00176 return pass;
00177 }
00178
00179 double FeketeQuadQuadrature::exact(int a, int b)
00180 {
00181 return 1.0 / (a + 1) / (b + 1);
00182 }
00183