Stokhos: Replace STEQR with PTEQR for issue trilinos#3542.

The LAPACK STEQR function, which computes eigenvalues/vectors of
symmetric, tridiagonal matrices, seg faults on the IBM Power systems for
some unknown reason.  However PTEQR, which does the same for SPD
systems, appears to work.  So I replaced to the call to STEQR in the
recurrence basis implementation with PTEQR.  This was more complicated
than it sounds because the recurrence matrix is not positive definite.
So I had to put in a shift to make it so for PTEQR.  This required
loosening a few tight numerical tolerances in some tests.

packages/stokhos/src/Stokhos_RecurrenceBasisImp.hpp

@@ -375,30 +375,59 @@ getQuadPoints(ordinal_type quad_order,
       normalizeRecurrenceCoefficients(a, b, c, d);
   }
 
-  // With normalized coefficients, A is symmetric and tri-diagonal, with
-  // diagonal = a, and off-diagonal = b, starting at b[1]
-  Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
-                                                                  num_points);
-  Teuchos::Array<value_type> workspace(2*num_points);
-  ordinal_type info_flag;
-  Teuchos::LAPACK<ordinal_type,value_type> my_lapack;
-
-  // compute the eigenvalues (stored in a) and right eigenvectors.
-  if (num_points == 1)
-    my_lapack.STEQR('I', num_points, &a[0], &b[0], eig_vectors.values(),
-                    num_points, &workspace[0], &info_flag);
-  else
-    my_lapack.STEQR('I', num_points, &a[0], &b[1], eig_vectors.values(),
-                    num_points, &workspace[0], &info_flag);
-
-  // eigenvalues are sorted by STEQR
   quad_points.resize(num_points);
   quad_weights.resize(num_points);
-  for (ordinal_type i=0; i<num_points; i++) {
-    quad_points[i] = a[i];
-    if (std::abs(quad_points[i]) < quad_zero_tol)
-      quad_points[i] = 0.0;
-    quad_weights[i] = beta[0]*eig_vectors[i][0]*eig_vectors[i][0];
+
+  if (num_points == 1) {
+    quad_points[0] = a[0];
+    quad_weights[0] = beta[0];
+  }
+  else {
+
+    // With normalized coefficients, A is symmetric and tri-diagonal, with
+    // diagonal = a, and off-diagonal = b, starting at b[1].  We use
+    // LAPACK'S PTEQR function to compute the eigenvalues/vectors of A
+    // to compute the quadrature points/weights respectively.  PTEQR requires
+    // an SPD matrix, so we need to shift the matrix to make it diagonally
+    // dominant (We could use STEQR which works for indefinite matrices, but
+    // this function has proven to be problematic on some platforms).
+    Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
+                                                                    num_points);
+    Teuchos::Array<value_type> workspace(4*num_points);
+    ordinal_type info_flag;
+    Teuchos::LAPACK<ordinal_type,value_type> my_lapack;
+
+    // Compute a shift to make matrix positive-definite
+    value_type max_a = 0.0;
+    value_type max_b = 0.0;
+    for (ordinal_type n = 0; n<num_points; n++) {
+      if (std::abs(a[n]) > max_a)
+        max_a = a[n];
+    }
+    for (ordinal_type n = 1; n<num_points; n++) {
+      if (std::abs(b[n]) > max_b)
+        max_b = b[n];
+    }
+    value_type shift = 1.0 + max_a + 2.0*max_b;
+
+    // Shift diagonal
+    for (ordinal_type n = 0; n<num_points; n++)
+      a[n] += shift;
+
+    // compute the eigenvalues (stored in a) and right eigenvectors.
+    my_lapack.PTEQR('I', num_points, &a[0], &b[1], eig_vectors.values(),
+                    num_points, &workspace[0], &info_flag);
+    TEUCHOS_TEST_FOR_EXCEPTION(info_flag != 0, std::logic_error,
+                               "PTEQR returned info = " << info_flag);
+
+    // (shifted) eigenvalues are sorted in descending order by PTEQR.
+    // Reorder them to ascending as in STEQR
+    for (ordinal_type i=0; i<num_points; i++) {
+      quad_points[i] = a[num_points-1-i]-shift;
+      if (std::abs(quad_points[i]) < quad_zero_tol)
+        quad_points[i] = 0.0;
+      quad_weights[i] = beta[0]*eig_vectors[num_points-1-i][0]*eig_vectors[num_points-1-i][0];
+    }
   }
 
   // Evalute basis at gauss points

packages/stokhos/test/UnitTest/Stokhos_HermiteBasisUnitTest.cpp

@@ -61,7 +61,7 @@ namespace HermiteBasisUnitTest {
     OrdinalType p;
     Stokhos::HermiteBasis<OrdinalType,ValueType> basis;
 
-    UnitTestSetup() : rtol(1e-12), atol(1e-7), p(10), basis(p) {}
+    UnitTestSetup() : rtol(1e-12), atol(1e-5), p(10), basis(p) {}
 
   };
 

packages/stokhos/test/UnitTest/Stokhos_LanczosUnitTest.cpp

@@ -89,7 +89,7 @@ struct Lanczos_PCE_Setup {
     func()
   {
     rtol = 1e-8;
-    atol = 1e-12;
+    atol = 1e-10;
     const OrdinalType d = 3;
     const OrdinalType p = 5;
 

packages/stokhos/test/UnitTest/Stokhos_NormalizedLegendreBasisUnitTest.cpp

@@ -61,7 +61,7 @@ namespace LegendreBasisUnitTest {
     OrdinalType p;
     Stokhos::LegendreBasis<OrdinalType,ValueType> basis;
 
-    UnitTestSetup() : rtol(1e-12), atol(1e-12), p(10), basis(p,true) {}
+    UnitTestSetup() : rtol(1e-12), atol(1e-10), p(10), basis(p,true) {}
 
   };
 

packages/stokhos/test/UnitTest/Stokhos_UnitTestHelpers.hpp

@@ -555,7 +555,7 @@ namespace Stokhos {
     out << "Discrete orthogonality error = " << x.normInf() << std::endl;
 
     // Compare PCE coefficients
-    bool success = Stokhos::compareSDM(x, "x", z, "zero", 1e-14, 1e-14, out);
+    bool success = Stokhos::compareSDM(x, "x", z, "zero", rel_tol, abs_tol, out);
 
     return success;
   }