Introduce cpp (and also adds lasr3 and steqr3)

.gitignore

@@ -43,3 +43,6 @@ build*
 DOCS/man
 DOCS/explore-html
 output_err
+
+# Editor config files
+.vscode/

LAPACK_CPP/Makefile

@@ -0,0 +1,45 @@
+CXX=g++
+# CXXFLAGS=-std=c++17 -Wall -Wextra -pedantic -g
+CXXFLAGS=-std=c++17 -O3 -march=native -ffast-math -g -DNDEBUG
+
+FC=gfortran
+FFLAGS=-Wall -Wextra -pedantic -g
+
+all: ./example/test_lasr3 ./example/profile_lasr3 ./example/profile_steqr3 $(OBJFILES)
+
+HEADERS = $(wildcard include/lapack_c/*.h) \
+		  $(wildcard include/lapack_c/*/*.h) \
+		  $(wildcard include/lapack_c/*/*/*.h) \
+		  $(wildcard include/lapack_cpp/*.hpp) \
+		  $(wildcard include/lapack_cpp/*/*.hpp) \
+		  $(wildcard include/lapack_cpp/*/*/*.hpp)
+
+SRCFILES = $(wildcard src/lapack_c/*.cpp) \
+		   $(wildcard src/lapack_c/*.f90) \
+		   $(wildcard src/lapack_cpp/*.cpp) \
+		   $(wildcard src/lapack_fortran/*.f90)
+
+OBJFILES = $(patsubst %.f90, %_f90.o ,$(patsubst %.cpp,%_cpp.o,$(SRCFILES)))
+
+# Fortran files
+src/lapack_fortran/%_f90.o: src/lapack_fortran/%.f90
+	$(FC) $(FFLAGS) -cpp -c -o $@ $<
+
+# C files
+
+src/lapack_c/%_cpp.o: src/lapack_c/%.cpp $(HEADERS)
+	$(CXX) $(CXXFLAGS) -I./include -c -o $@ $<
+
+src/lapack_c/%_f90.o: src/lapack_c/%.f90
+	$(FC) $(FFLAGS) -cpp -c -o $@ $<
+
+# C++ files
+src/lapack_cpp/%_cpp.o: src/lapack_cpp/%.cpp $(HEADERS)
+	$(CXX) $(CXXFLAGS) -I./include -c -o $@ $<
+
+# Test files
+example/%: example/%.cpp $(OBJFILES) $(HEADERS)
+	$(CXX) $(CXXFLAGS) -o $@ $< $(OBJFILES) -I./include -lstdc++ -lgfortran -lblas -llapack
+
+clean:
+	find . -type f -name '*.o' -delete

LAPACK_CPP/example/.gitignore

@@ -0,0 +1,3 @@
+test_lasr3
+profile_lasr3
+profile_steqr3

LAPACK_CPP/example/profile_lasr3.cpp

@@ -0,0 +1,111 @@
+
+#include <chrono> // for high_resolution_clock
+#include <complex>
+#include <iostream>
+
+#include "lapack_cpp.hpp"
+
+using namespace lapack_cpp;
+
+template <typename T,
+          Layout layout = Layout::ColMajor,
+          typename idx_t = lapack_idx_t>
+void profile_lasr3(idx_t m, idx_t n, idx_t k, Side side, Direction direction)
+{
+    idx_t n_rot = side == Side::Left ? m - 1 : n - 1;
+
+    MemoryBlock<T, idx_t> A_(m, n, layout);
+    Matrix<T, layout, idx_t> A(m, n, A_);
+
+    MemoryBlock<real_t<T>, idx_t> C_(n_rot, k, layout);
+    Matrix<real_t<T>, layout, idx_t> C(n_rot, k, C_);
+
+    MemoryBlock<T, idx_t> S_(n_rot, k, layout);
+    Matrix<T, layout, idx_t> S(n_rot, k, S_);
+
+    randomize(A);
+
+    // Generate random, but valid rotations
+    randomize(C);
+    randomize(S);
+    for (idx_t i = 0; i < n_rot; ++i)
+    {
+        for (idx_t j = 0; j < k; ++j)
+        {
+            T f = C(i, j);
+            T g = S(i, j);
+            T r;
+            lartg(f, g, C(i, j), S(i, j), r);
+        }
+    }
+
+    MemoryBlock<T, idx_t> A_copy_(m, n, layout);
+    Matrix<T, layout, idx_t> A_copy(m, n, A_copy_);
+
+    A_copy = A;
+
+    MemoryBlock<T, idx_t> work(
+        lasr3_workquery(side, direction, C.as_const(), S.as_const(), A));
+
+    const idx_t n_timings = 100;
+    const idx_t n_warmup = 20;
+    std::vector<float> timings(n_timings);
+
+    for (idx_t i = 0; i < n_timings; ++i)
+    {
+        A = A_copy;
+        auto start = std::chrono::high_resolution_clock::now();
+        lasr3(side, direction, C.as_const(), S.as_const(), A, work);
+        auto end = std::chrono::high_resolution_clock::now();
+        timings[i] = std::chrono::duration<float>(end - start).count();
+    }
+
+    float mean = 0;
+    for (idx_t i = n_warmup; i < n_timings; ++i)
+        mean += timings[i];
+    mean /= n_timings - n_warmup;
+
+    float std_dev = 0;
+    for (idx_t i = n_warmup; i < n_timings; ++i)
+        std_dev += (timings[i] - mean) * (timings[i] - mean);
+    std_dev = std::sqrt(std_dev / (n_timings - n_warmup - 1));
+
+    long nflops =
+        side == Side::Left ? 6. * (m - 1) * n * k : 6. * (n - 1) * m * k;
+
+    std::cout << "m = " << m << ", n = " << n << ", k = " << k
+              << ", side = " << (char)side
+              << ", direction = " << (char)direction << ", mean time = " << mean
+              << " s"
+              << ", std dev = " << std_dev / mean * 100 << " %"
+              << ", flop rate = " << nflops / mean * 1.0e-9 << " GFlops"
+              << std::endl;
+}
+
+int main()
+{
+    typedef lapack_idx_t idx_t;
+    typedef double T;
+
+    const idx_t nb = 1000;
+    const idx_t k = 64;
+
+    for (int s = 0; s < 2; ++s)
+    {
+        Side side = s == 0 ? Side::Left : Side::Right;
+        for (int d = 0; d < 2; ++d)
+        {
+            Direction direction =
+                d == 0 ? Direction::Forward : Direction::Backward;
+            for (idx_t n_rot = 99; n_rot <= 1600; n_rot += 100)
+            {
+                idx_t m = side == Side::Left ? n_rot + 1 : nb;
+                idx_t n = side == Side::Left ? nb : n_rot + 1;
+
+                profile_lasr3<T, Layout::ColMajor, lapack_idx_t>(m, n, k, side,
+                                                                 direction);
+            }
+        }
+    }
+    return 0;
+}

LAPACK_CPP/example/profile_steqr3.cpp

@@ -0,0 +1,102 @@
+
+#include <chrono> // for high_resolution_clock
+#include <complex>
+#include <iostream>
+
+#include "lapack_cpp.hpp"
+
+using namespace lapack_cpp;
+
+template <typename T,
+          Layout layout = Layout::ColMajor,
+          typename idx_t = lapack_idx_t>
+void profile_steqr3(idx_t n, bool use_fortran = false)
+{
+    MemoryBlock<T, idx_t> Z_(n, n, layout);
+    Matrix<T, layout, idx_t> Z(n, n, Z_);
+
+    MemoryBlock<T, idx_t> d_(n);
+    Vector<T, idx_t> d(n, d_);
+
+    MemoryBlock<T, idx_t> e_(n - 1);
+    Vector<T, idx_t> e(n - 1, e_);
+
+    MemoryBlock<T, idx_t> d_copy_(n);
+    Vector<T, idx_t> d_copy(n, d_copy_);
+
+    MemoryBlock<T, idx_t> e_copy_(n - 1);
+    Vector<T, idx_t> e_copy(n - 1, e_copy_);
+
+    randomize(d);
+    randomize(e);
+
+    d_copy = d;
+    e_copy = e;
+
+    CompQ compz = CompQ::Initialize;
+
+    const idx_t n_timings = 100;
+    const idx_t n_warmup = 50;
+    std::vector<float> timings(n_timings);
+
+    if (use_fortran)
+    {
+        MemoryBlock<real_t<T>, idx_t> rwork(2 * n - 2);
+        for (idx_t i = 0; i < n_timings; ++i)
+        {
+            d = d_copy;
+            e = e_copy;
+            auto start = std::chrono::high_resolution_clock::now();
+            steqr(compz, d, e, Z, rwork);
+            auto end = std::chrono::high_resolution_clock::now();
+            timings[i] = std::chrono::duration<float>(end - start).count();
+            std::cout << i << " " << timings[i] << '\r' << std::flush;
+        }
+    }
+    else
+    {
+        MemoryBlock<T, idx_t> work(steqr3_workquery(compz, d, e, Z));
+        MemoryBlock<real_t<T>, idx_t> rwork(steqr3_rworkquery(compz, d, e, Z));
+        for (idx_t i = 0; i < n_timings; ++i)
+        {
+            d = d_copy;
+            e = e_copy;
+            auto start = std::chrono::high_resolution_clock::now();
+            steqr3(compz, d, e, Z, work, rwork);
+            auto end = std::chrono::high_resolution_clock::now();
+            timings[i] = std::chrono::duration<float>(end - start).count();
+            std::cout << i << " " << timings[i] << '\r' << std::flush;
+        }
+    }
+
+    float mean = 0;
+    for (idx_t i = n_warmup; i < n_timings; ++i)
+        mean += timings[i];
+    mean /= n_timings - n_warmup;
+
+    float std_dev = 0;
+    for (idx_t i = n_warmup; i < n_timings; ++i)
+        std_dev += (timings[i] - mean) * (timings[i] - mean);
+    std_dev = std::sqrt(std_dev / (n_timings - n_warmup - 1));
+
+    std::cout << "n = " << n << ", using fortran: " << use_fortran << ", mean time = " << mean
+              << " s"
+              << ", std dev = " << std_dev / mean * 100 << " %"
+              << std::endl;
+}
+
+int main()
+{
+    typedef lapack_idx_t idx_t;
+    typedef double T;
+
+    for (int n = 32; n <= 2000; n *= 2)
+    {
+        profile_steqr3<T, Layout::ColMajor, idx_t>(n, false);
+    }
+    for (int n = 32; n <= 2000; n *= 2)
+    {
+        profile_steqr3<T, Layout::ColMajor, idx_t>(n, true);
+    }
+    return 0;
+}

LAPACK_CPP/example/test_lasr3.cpp

@@ -0,0 +1,115 @@
+
+#include <complex>
+#include <iostream>
+
+#include "lapack_cpp.hpp"
+
+using namespace lapack_cpp;
+
+template <typename T,
+          Layout layout = Layout::ColMajor,
+          typename idx_t = lapack_idx_t>
+void test_lasr3(idx_t m, idx_t n, idx_t k, Side side, Direction direction)
+{
+    idx_t n_rot = side == Side::Left ? m - 1 : n - 1;
+
+    MemoryBlock<T, idx_t> A_(m, n, layout);
+    Matrix<T, layout, idx_t> A(m, n, A_);
+
+    MemoryBlock<real_t<T>, idx_t> C_(n_rot, k, layout);
+    Matrix<real_t<T>, layout, idx_t> C(n_rot, k, C_);
+
+    MemoryBlock<T, idx_t> S_(n_rot, k, layout);
+    Matrix<T, layout, idx_t> S(n_rot, k, S_);
+
+    randomize(A);
+
+    // Generate random, but valid rotations
+    randomize(C);
+    randomize(S);
+    for (idx_t i = 0; i < n_rot; ++i) {
+        for (idx_t j = 0; j < k; ++j) {
+            T f = C(i, j);
+            T g = S(i, j);
+            T r;
+            lartg(f, g, C(i, j), S(i, j), r);
+        }
+    }
+
+    MemoryBlock<T, idx_t> A_copy_(m, n, layout);
+    Matrix<T, layout, idx_t> A_copy(m, n, A_copy_);
+
+    A_copy = A;
+
+    // Apply rotations using simple loops as test
+    if (side == Side::Left) {
+        if (direction == Direction::Forward) {
+            for (idx_t i = 0; i < k; ++i)
+                for (idx_t j = 0; j < n_rot; ++j)
+                    rot(A_copy.row(j), A_copy.row(j + 1), C(j, i), S(j, i));
+        }
+        else {
+            for (idx_t i = 0; i < k; ++i)
+                for (idx_t j = n_rot - 1; j >= 0; --j)
+                    rot(A_copy.row(j), A_copy.row(j + 1), C(j, i), S(j, i));
+        }
+    }
+    else {
+        if (direction == Direction::Forward) {
+            for (idx_t i = 0; i < k; ++i)
+                for (idx_t j = 0; j < n_rot; ++j)
+                    rot(A_copy.column(j), A_copy.column(j + 1), C(j, i),
+                        conj(S(j, i)));
+        }
+        else {
+            for (idx_t i = 0; i < k; ++i)
+                for (idx_t j = n_rot - 1; j >= 0; --j)
+                    rot(A_copy.column(j), A_copy.column(j + 1), C(j, i),
+                        conj(S(j, i)));
+        }
+    }
+
+    lasr3(side, direction, C.as_const(), S.as_const(), A);
+
+    // Check that the result is the same
+    for (idx_t i = 0; i < m; ++i)
+        for (idx_t j = 0; j < n; ++j)
+            A_copy(i, j) -= A(i, j);
+
+    real_t<T> err = 0.;
+    for (idx_t i = 0; i < m; ++i)
+        for (idx_t j = 0; j < n; ++j)
+            err = std::max(err, abs(A_copy(i, j)));
+
+    // This test should really be relative
+    if( err > 1e-5 ){
+        std::cout << "Failed test_lasr3 with parameters: " << m << ", " << n << ", " << k << ", " << (char)side << ", " << (char)direction << std::endl;
+        std::cout << "Error: " << err << std::endl;
+    } else {
+        std::cout << "Passed test_lasr3 with parameters: " << m << ", " << n << ", " << k << ", " << (char)side << ", " << (char)direction << std::endl;
+    }
+
+    // print(A_copy);
+}
+
+int main()
+{
+    typedef lapack_idx_t idx_t;
+    typedef std::complex<float> T;
+
+    for (idx_t nb = 64; nb <= 2000; nb *= 2) {
+        for (idx_t n_rot = 2; n_rot <= 2000; n_rot*=2) {
+            for (idx_t k = 32; k <= 128; k+=32) {
+                test_lasr3<T, Layout::ColMajor, lapack_idx_t>(
+                    n_rot + 1, nb, k, Side::Left, Direction::Forward);
+                test_lasr3<T, Layout::ColMajor, lapack_idx_t>(
+                    nb, n_rot + 1, k, Side::Right, Direction::Forward);
+                test_lasr3<T, Layout::ColMajor, lapack_idx_t>(
+                    n_rot + 1, nb, k, Side::Left, Direction::Backward);
+                test_lasr3<T, Layout::ColMajor, lapack_idx_t>(
+                    nb, n_rot + 1, k, Side::Right, Direction::Backward);
+            }
+        }
+    }
+    return 0;
+}