Add linear algebra file

CMakeLists.txt

@@ -7,6 +7,8 @@ include_directories(${CPLEX_INCLUDE_DIR})
 
 find_package(GMP REQUIRED)
 find_package(OpenMP REQUIRED)
+find_package(BLAS REQUIRED)
+find_package(LAPACKE REQUIRED)
 
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Werror ${OpenMP_CXX_FLAGS} -O3")

cmake/FindLAPACKE.cmake

@@ -0,0 +1,11 @@
+find_library(LAPACKE_LIBRARIES
+        NAMES lapacke)
+
+find_path(LAPACKE_INCLUDE_DIR NAMES lapacke.h)
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(LAPACKE
+        DEFAULT_MSG
+        LAPACKE_LIBRARIES
+        LAPACKE_INCLUDE_DIR)
+mark_as_advanced(LAPACKE_LIBRARIES LAPACKE_INCLUDE_DIR)

multirow/CMakeLists.txt

@@ -1,33 +1,43 @@
 set(COMMON_SOURCES
-    src/cg.c
+        src/cg.c
-    src/double.c
+        src/double.c
-    src/geometry.c
+        src/geometry.c
-    src/lfree2d.c
+        src/lfree2d.c
-    src/lp.c
+        src/lp.c
-    src/mir.c
+        src/mir.c
-    src/util.c
+        src/util.c
-    src/rational.c
+        src/rational.c
-    src/stats.c
+        src/stats.c
-    include/multirow/cg.h
+        src/linalg.c
-    include/multirow/double.h
+        include/multirow/cg.h
-    include/multirow/geometry.h
+        include/multirow/double.h
-    include/multirow/lfree2d.h
+        include/multirow/geometry.h
-    include/multirow/lp.h
+        include/multirow/lfree2d.h
-    include/multirow/mir.h
+        include/multirow/lp.h
-    include/multirow/rational.h
+        include/multirow/mir.h
-    include/multirow/stats.h
+        include/multirow/rational.h
-    include/multirow/params.h
+        include/multirow/stats.h
-    include/multirow/util.h)
+        include/multirow/params.h
+        include/multirow/util.h
+        include/multirow/linalg.h)
 
 set(TEST_SOURCES
-    tests/double-test.cpp
+        tests/double-test.cpp
-    tests/cg-test.cpp
+        tests/cg-test.cpp
-    tests/geometry-test.cpp)
+        tests/geometry-test.cpp
+        tests/linalg-test.cpp)
 
 add_library(multirow_static ${COMMON_SOURCES})
 set_target_properties(multirow_static PROPERTIES OUTPUT_NAME lifting)
-target_link_libraries(multirow_static ${CPLEX_LIBRARIES})
+target_link_libraries(multirow_static
-target_include_directories (multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+        ${CPLEX_LIBRARIES}
+        ${BLAS_LIBRARIES}
+        ${LAPACKE_LIBRARIES})
+target_include_directories(multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 add_executable(multirow-test.run ${COMMON_SOURCES} ${TEST_SOURCES})
-target_link_libraries(multirow-test.run ${CPLEX_LIBRARIES} gtest_main)
+target_link_libraries(multirow-test.run
+        ${CPLEX_LIBRARIES}
+        ${BLAS_LIBRARIES}
+        ${LAPACKE_LIBRARIES}
+        gtest_main)

multirow/include/multirow/linalg.h

@@ -0,0 +1,40 @@
+/* Copyright (c) 2015-2017 Alinson Xavier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+/**
+ * Receives an n-dimensional vector x, a scalar alpha and sets x <- alpha * x
+ */
+void LINALG_scale(int n, double *x, double alpha);
+
+/**
+ * Receives two n-dimensional vectors x and y and returns the dot product of x
+ * and y.
+ */
+double LINALG_dot(int n, double *x, double *y);
+
+/**
+ * Receives an n-dimensional vector x and returns the 1-norm of x.
+ */
+double LINALG_norm(int n, double *x);
+
+/**
+ * Given an invertible n-by-n matrix A and an n-dimensional vector b, this
+ * function set x to A^(-1) b. Returns zero if the operation is successful and
+ * non-zero otherwise.
+ */
+int LINALG_solve(int n, double *A, double *b, double *x);

multirow/src/linalg.c

@@ -0,0 +1,53 @@
+/*
+ * Copyright (C) 2016 Álinson Santos Xavier <isoron@gmail.com>
+ *
+ * This file is part of Loop Habit Tracker.
+ *
+ * Loop Habit Tracker is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by the
+ * Free Software Foundation, either version 3 of the License, or (at your
+ * option) any later version.
+ *
+ * Loop Habit Tracker is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+ * more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include<cblas.h>
+#include<lapacke.h>
+#include <multirow/util.h>
+
+void LINALG_scale(int n, double *x, double alpha)
+{
+    cblas_dscal(n, alpha, x, 1);
+}
+
+double LINALG_dot(int n, double *x, double *y)
+{
+    return cblas_ddot(n, x, 1, y, 1);
+}
+
+double LINALG_norm(int n, double *x)
+{
+    return cblas_dasum(n, x, 1);
+}
+
+int LINALG_solve(int n, double *A, double *b, double *x)
+{
+    int rval = 0;
+    int ignored[n];
+    double A_copy[n * n];
+
+    memcpy(x, b, n * sizeof(double));
+    memcpy(A_copy, A, n * n * sizeof(double));
+
+    rval = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A_copy, n, ignored, x, 1);
+    abort_if(rval, "LAPACKE_dgesv failed");
+
+CLEANUP:
+    return rval;
+}

multirow/tests/linalg-test.cpp

@@ -0,0 +1,83 @@
+/* Copyright (c) 2015-2017 Alinson Xavier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gtest/gtest.h>
+
+extern "C" {
+#include <multirow/util.h>
+#include <multirow/linalg.h>
+}
+
+const double E = 1e-6;
+
+TEST(LinAlgTest, scale_test)
+{
+    double x[] = {1.0, 2.0, 3.0};
+    LINALG_scale(3, x, 2.0);
+    EXPECT_NEAR(x[0], 2.0, E);
+    EXPECT_NEAR(x[1], 4.0, E);
+    EXPECT_NEAR(x[2], 6.0, E);
+}
+
+TEST(LinAlgTest, dot_test)
+{
+    double x[] = { 1.0, 2.0, 3.0 };
+    double y[] = { 3.0, 4.0, 5.0 };
+    double dot = LINALG_dot(3, x, y);
+    EXPECT_NEAR(dot, 26.0, E);
+}
+
+TEST(LinAlgTest, norm_test)
+{
+    double x[] = { 1.0,  2.0, -3.0 };
+    double y[] = { 3.0, -4.0, -5.0 };
+    double x_norm = LINALG_norm(3, x);
+    double y_norm = LINALG_norm(3, y);
+    EXPECT_NEAR(x_norm, 6.0, E);
+    EXPECT_NEAR(y_norm, 12.0, E);
+}
+
+TEST(LinAlgTest, solve_test)
+{
+    int rval = 0;
+
+    double A[] = {
+            2.0, 1.0,  3.0,
+            2.0, 6.0,  8.0,
+            6.0, 8.0, 18.0,
+    };
+    double b[] = { 1.0, 3.0, 5.0 };
+    double x[] = { 0.0, 0.0, 0.0 };
+
+    rval = LINALG_solve(3, A, b, x);
+    abort_if(rval, "LINALG_solve failed");
+
+    // Should compute x correctly
+    EXPECT_NEAR(x[0], 0.3, E);
+    EXPECT_NEAR(x[1], 0.4, E);
+    EXPECT_NEAR(x[2], 0.0, E);
+
+    // Should not modify A and b
+    EXPECT_EQ(A[0], 2.0);
+    EXPECT_EQ(A[1], 1.0);
+    EXPECT_EQ(A[2], 3.0);
+    EXPECT_EQ(b[0], 1.0);
+    EXPECT_EQ(b[1], 3.0);
+    EXPECT_EQ(b[2], 5.0);
+
+CLEANUP:
+    if(rval) FAIL();
+}