345 Teuchos::Array<value_type>& quad_points,
346 Teuchos::Array<value_type>& quad_weights,
347 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const
354 ordinal_type num_points =
355 static_cast<ordinal_type
>(std::ceil((quad_order+1)/2.0));
356 Teuchos::Array<value_type> a(num_points,0);
357 Teuchos::Array<value_type> b(num_points,0);
358 Teuchos::Array<value_type> c(num_points,0);
359 Teuchos::Array<value_type> d(num_points,0);
362 if(num_points > p+1){
363 bool is_normalized = computeRecurrenceCoefficients(num_points, a, b, c, d);
365 normalizeRecurrenceCoefficients(a, b, c, d);
368 for(ordinal_type n = 0; n<num_points; n++){
375 normalizeRecurrenceCoefficients(a, b, c, d);
378 quad_points.resize(num_points);
379 quad_weights.resize(num_points);
381 if (num_points == 1) {
382 quad_points[0] = a[0];
383 quad_weights[0] = beta[0];
394 Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
396 Teuchos::Array<value_type> workspace(4*num_points);
397 ordinal_type info_flag;
401 value_type max_a = 0.0;
402 value_type max_b = 0.0;
403 for (ordinal_type n = 0; n<num_points; n++) {
404 if (std::abs(a[n]) > max_a)
407 for (ordinal_type n = 1; n<num_points; n++) {
408 if (std::abs(b[n]) > max_b)
411 value_type shift = 1.0 + max_a + 2.0*max_b;
414 for (ordinal_type n = 0; n<num_points; n++)
418 my_lapack.PTEQR(
'I', num_points, &a[0], &b[1], eig_vectors.values(),
419 num_points, &workspace[0], &info_flag);
420 TEUCHOS_TEST_FOR_EXCEPTION(info_flag != 0, std::logic_error,
421 "PTEQR returned info = " << info_flag);
425 for (ordinal_type i=0; i<num_points; i++) {
426 quad_points[i] = a[num_points-1-i]-shift;
427 if (std::abs(quad_points[i]) < quad_zero_tol)
428 quad_points[i] = 0.0;
429 quad_weights[i] = beta[0]*eig_vectors[num_points-1-i][0]*eig_vectors[num_points-1-i][0];
434 quad_values.resize(num_points);
435 for (ordinal_type i=0; i<num_points; i++) {
436 quad_values[i].resize(p+1);
437 evaluateBases(quad_points[i], quad_values[i]);
495 Teuchos::Array<value_type>& a,
496 Teuchos::Array<value_type>& b,
497 Teuchos::Array<value_type>& c,
498 Teuchos::Array<value_type>&
g)
const
500 ordinal_type n = a.size();
502 b[0] = std::sqrt(b[0]);
504 for (ordinal_type k=1; k<n; k++) {
506 b[k] = std::sqrt((b[k]*
g[k])/(c[k]*c[k-1]));
509 for (ordinal_type k=0; k<n; k++)
510 c[k] = value_type(1);