Test: add xUNBDB5 test for LAPACK issue Reference-LAPACK#634

TESTING/cpp/CMakeLists.txt

@@ -47,6 +47,7 @@ add_cxx_test(tools_tests tools_tests.cpp)
 add_cxx_test(xGGQRCS_tests xGGQRCS_tests.cpp)
 add_cxx_test(xLASRTI_tests xLASRTI_tests.cpp)
 add_cxx_test(xLASRTR_tests xLASRTR_tests.cpp)
+add_cxx_test(xUNBDB5_tests xUNBDB5_tests.cpp)
 add_cxx_test(xUNCSD2BY1_tests xUNCSD2BY1_tests.cpp)
 
 

TESTING/cpp/xUNBDB5_tests.cpp

@@ -0,0 +1,88 @@
+/*
+ * Copyright (c) 2021 Christoph Conrads
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *  * Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ *  * Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *  * Neither the name of the copyright holders nor the
+ *    names of its contributors may be used to endorse or promote products
+ *    derived from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "config.hpp"
+#include "tools.hpp"
+
+#include <boost/assert.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/test/unit_test.hpp>
+#include <cstddef>
+
+using Integer = lapack::integer_t;
+
+namespace ublas = boost::numeric::ublas;
+namespace tools = lapack::tools;
+
+template<typename Number,
+         typename std::enable_if<std::is_fundamental<Number>::value, int>::type
+             * = nullptr>
+void miscomputed_u1_impl(Number) {
+    using Real = typename tools::real_from<Number>::type;
+    using Matrix = ublas::matrix<Number, ublas::column_major>;
+
+    auto m = std::size_t{2};
+    auto n = std::size_t{3};
+    auto p = std::size_t{2};
+    auto ldq = m + p;
+    auto Q = Matrix(ldq, n);
+    constexpr auto d = std::numeric_limits<Real>::digits - 1;
+    constexpr auto eps = std::numeric_limits<Real>::epsilon();
+    // the exact value does matter, it should just be signficantly smaller than
+    // the machine epsilon the preclude round-off errors
+    auto mu = std::ldexp(eps, -d);
+    // clang-format off
+    Q(0, 0) = 0;  Q(0, 1) = 1;  Q(0, 2) = -mu;
+    Q(1, 0) = 0;  Q(1, 1) = 0;  Q(1, 2) = +mu;
+    Q(2, 0) = 1;  Q(2, 1) = 0;  Q(2, 2) = 0;
+    Q(3, 0) = 0;  Q(3, 1) = mu; Q(3, 2) = 1;
+    // clang-format on
+
+    BOOST_VERIFY(tools::is_almost_isometric(Q));
+
+    auto nan = tools::not_a_number<Number>::value;
+    auto x = ublas::vector<Number>(m + p, Number{});
+    auto lwork = n;
+    auto work = ublas::vector<Number>(lwork, nan);
+    auto ret = lapack::xUNBDB5(m, p, n, &x(0), 1, &x(m), 1, &Q(0, 0), ldq,
+                               &Q(m, 0), ldq, &work(0), lwork);
+
+    auto norm_x = lapack::xNRM2(x.size(), &x(0), 1);
+    BOOST_REQUIRE_EQUAL(ret, 0);
+    // for certain values of mu, ublas::norm_2 returns zero, e.g., for
+    // eps*2^-20 in single precision
+    BOOST_REQUIRE_GT(norm_x, 0);
+    BOOST_CHECK_LE(std::abs(x(1) - norm_x), 2 * norm_x * eps);
+    std::printf("%16.10e %16.10e %16.10e %16.10e\n", x(0), x(1), x(2), x(3));
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(miscomputed_u1, Number,
+                              lapack::supported_real_types) {
+    miscomputed_u1_impl(Number{});
+}