Introduce cpp (and also adds lasr3 and steqr3)

LAPACK_CPP/Makefile

@@ -0,0 +1,51 @@
+CXX=g++
+# CXXFLAGS=-std=c++17 -Wall -Wextra -pedantic -g -Og
+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
+
+HEADERS = $(wildcard lapack_c/include/*.h) \
+		  $(wildcard lapack_c/include/*/*.h) \
+		  $(wildcard lapack_c/include/*/*/*.h) \
+		  $(wildcard lapack_cpp/include/*.hpp) \
+		  $(wildcard lapack_cpp/include/*/*.hpp) \
+		  $(wildcard lapack_cpp/include/*/*/*.hpp)
+
+OBJFILES = src/lapack_cpp/rot_cpp.o \
+		   src/lapack_cpp/lartg_cpp.o \
+		   src/lapack_c/rot_c.o \
+		   src/lapack_c/lartg_c.o \
+		   src/lapack_c/lasr3_c.o \
+		   src/lapack_fortran/slasr3.o \
+		   src/lapack_fortran/dlasr3.o \
+		   src/lapack_fortran/clasr3.o \
+		   src/lapack_fortran/zlasr3.o \
+
+# Fortran files
+src/lapack_fortran/%.o: src/lapack_fortran/%.f90
+	$(FC) $(FFLAGS) -cpp -c -o $@ $<
+
+# C files
+
+src/lapack_c/%.o: src/lapack_c/%.cpp
+	$(CXX) $(CXXFLAGS) -I./include -c -o $@ $<
+
+src/lapack_c/%.o: src/lapack_c/%.f90
+	$(FC) $(FFLAGS) -c -o $@ $<
+
+# C++ files
+src/lapack_cpp/%.o: src/lapack_cpp/%.cpp $(HEADERS)
+	$(CXX) $(CXXFLAGS) -I./include -c -o $@ $<
+
+# Test files
+example/test_lasr3: example/test_lasr3.cpp $(OBJFILES) $(HEADERS)
+	$(CXX) $(CXXFLAGS) -o $@ $< $(OBJFILES) -I./include -lstdc++ -lgfortran -lblas -llapack
+
+example/profile_lasr3: example/profile_lasr3.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,2 @@
+test_lasr3
+profile_lasr3

LAPACK_CPP/example/profile_lasr3.cpp

@@ -0,0 +1,99 @@
+
+#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;
+
+    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"
+              << ", 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/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;
+}

LAPACK_CPP/include/lapack_c.h

@@ -0,0 +1,13 @@
+#ifndef LAPACK_C_HPP
+#define LAPACK_C_HPP
+
+#include "lapack_c/util.h"
+
+// BLAS
+#include "lapack_c/rot.h"
+
+// LAPACK
+#include "lapack_c/lartg.h"
+
+
+#endif

LAPACK_CPP/include/lapack_c/lartg.h

@@ -0,0 +1,22 @@
+#ifndef LAPACK_C_LARTG_HPP
+#define LAPACK_C_LARTG_HPP
+
+#include "lapack_c/util.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif  // __cplusplus
+
+void lapack_c_slartg(float f, float g, float* cs, float* sn, float* r);
+
+void lapack_c_dlartg(double f, double g, double* cs, double* sn, double* r);
+
+void lapack_c_clartg(lapack_float_complex f, lapack_float_complex g, float* cs, lapack_float_complex* sn, lapack_float_complex* r);
+
+void lapack_c_zlartg(lapack_double_complex f, lapack_double_complex g, double* cs, lapack_double_complex* sn, lapack_double_complex* r);
+
+#ifdef __cplusplus
+}
+#endif  // __cplusplus
+
+#endif