diff --git a/CMakeLists.txt b/CMakeLists.txt index 09908e4..043c180 100644 --- a/CMakeLists.txt +++ b/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") diff --git a/cmake/FindLAPACKE.cmake b/cmake/FindLAPACKE.cmake new file mode 100644 index 0000000..4e827cd --- /dev/null +++ b/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) diff --git a/multirow/CMakeLists.txt b/multirow/CMakeLists.txt index ff7a311..d7314ea 100644 --- a/multirow/CMakeLists.txt +++ b/multirow/CMakeLists.txt @@ -1,33 +1,43 @@ set(COMMON_SOURCES - src/cg.c - src/double.c - src/geometry.c - src/lfree2d.c - src/lp.c - src/mir.c - src/util.c - src/rational.c - src/stats.c - include/multirow/cg.h - include/multirow/double.h - include/multirow/geometry.h - include/multirow/lfree2d.h - include/multirow/lp.h - include/multirow/mir.h - include/multirow/rational.h - include/multirow/stats.h - include/multirow/params.h - include/multirow/util.h) + src/cg.c + src/double.c + src/geometry.c + src/lfree2d.c + src/lp.c + src/mir.c + src/util.c + src/rational.c + src/stats.c + src/linalg.c + include/multirow/cg.h + include/multirow/double.h + include/multirow/geometry.h + include/multirow/lfree2d.h + include/multirow/lp.h + include/multirow/mir.h + include/multirow/rational.h + include/multirow/stats.h + include/multirow/params.h + include/multirow/util.h + include/multirow/linalg.h) set(TEST_SOURCES - tests/double-test.cpp - tests/cg-test.cpp - tests/geometry-test.cpp) + tests/double-test.cpp + tests/cg-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_include_directories (multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(multirow_static + ${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) diff --git a/multirow/include/multirow/linalg.h b/multirow/include/multirow/linalg.h new file mode 100644 index 0000000..142aca8 --- /dev/null +++ b/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 . + */ + +#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); \ No newline at end of file diff --git a/multirow/src/linalg.c b/multirow/src/linalg.c new file mode 100644 index 0000000..21b34ee --- /dev/null +++ b/multirow/src/linalg.c @@ -0,0 +1,53 @@ +/* + * Copyright (C) 2016 Álinson Santos Xavier + * + * 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 . + */ + +#include +#include +#include + +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; +} \ No newline at end of file diff --git a/multirow/tests/linalg-test.cpp b/multirow/tests/linalg-test.cpp new file mode 100644 index 0000000..2ddafec --- /dev/null +++ b/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 . + */ + +#include + +extern "C" { +#include +#include +} + +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(); +} \ No newline at end of file