Commit 131ce181 by Pratyuksh Bansal

### Merge branch 'master' of gitlab.math.ethz.ch:NumCSE/NumCSE

parents ee305b1e 999e80aa
 #include #include #include using namespace Eigen; /* @brief Solve the system Ry=c * for the upper triangular matrix R * This could help you in your implementation * of solve_LSE() * \param[in] R nxn regular, upper triangular matrix * \param[in] c n dim vector * \return y n dim result vector */ /* SAM_LISTING_BEGIN_2 */ VectorXd solve_R(const MatrixXd& R, const VectorXd& c) { int n = R.rows(); assert(n == R.cols() && n == c.size() && "Input dimensions must agree"); // Initialize VectorXd y(n); #if SOLUTION // Since R is upper triangular, we can solve by backwards substitution for (int i = n-1; i >= 0; --i) { y(i) = c(i); for (int j = n-1; j > i ; --j) { y(i) -= R(i,j) * y(j); } y(i) /= R(i,i); } #else // Implementing this function could help you in solve_LSE() #endif return y; } /* SAM_LISTING_END_2 */ /* @brief Solve the System Ax=b * for A << R, v, * u.transpose(), 0; * \param[in] R nxn regular, upper triangular matrix * \param[in] v n dim vector * \param[in] u n dim vector * \param[in] b n+1 dim vector * \return x n+1 dim result vector */ /* SAM_LISTING_BEGIN_1 */ VectorXd solve_LSE(const MatrixXd& R, const VectorXd& v, const VectorXd& u, const VectorXd& b) { unsigned n = R.rows(); assert(R.cols() == n && "R has to be square"); assert(n == v.size() && n == u.size() && n+1 == b.size() && "Input dimensions must agree"); // Initialize VectorXd y(n+1), x(n+1); #if SOLUTION // Solve the system Ax=b by LU-Decomposition. // Solve Ly = b through forward substitution. // Due to the special structure of our L, // the first n entries of y are easy: y.head(n) = b.head(n); // The last element of y is given by $y_n = b_n - u^T\mathbf{R}^{-1}y_{0...n-1}$ y(n) = b(n) - u.transpose() * solve_R(R, y.head(n)); // Solve Ux = y by backward substitution. // First we build U MatrixXd U(n+1,n+1); U << R, v, VectorXd::Zero(n).transpose(), -u.transpose()*solve_R(R,v); // Then we solve Ux = y x = solve_R(U,y); // \iffalse Latex comment-delimiter so this part doesn't appear // in the solution and screws up formatting // Note this could be done using less memory by not constructing // U, doing the first step of the back substitution "by hand" and then // calling solve_R() // VectorXd x(n+1); // x(n) = y(n) / (-u.transpose()*solve_R(R,v)); // x.head(n) = solve_R(R,y.head(n) - v*x(n)); //\fi #else // Solve the LSE using LU-decomposition and the expression // for L and U that you derived #endif return x; } /* SAM_LISTING_END_1 */ int main() { // Vectors for testing unsigned n = 10; VectorXd v,u,b; u = v = VectorXd::Random(n); b = VectorXd::Random(n+1); // Upper triangular matrix MatrixXd R(n,n); for (unsigned i = 0; i < n; ++i) { for (unsigned j = i; j < n; ++j) { R(i,j) = rand(); //Bad RNG, but sufficient here } } R /= RAND_MAX; //"norm" R for numerical stability // Build matrix A for Eigensolver MatrixXd A(n+1,n+1); A << R, v, u.transpose(), 0; double error = (solve_LSE(R,v,u,b) - A.colPivHouseholderQr().solve(b)).norm(); if (error > 1e-8) { std::cout << "solve_LSE() returns a different solution than Eigen" << std::endl; } else { std::cout << "solve_LSE() and Eigen get the same result" << std::endl; } }
 { "author" : "ochsnerd", "email" : "", "contributors" : "", "name": "BlockLU", "all" : ["blockLU.cpp"] }
 project(BlockLU) cmake_minimum_required(VERSION 2.8) # # setup compiler # # use c++11 if(CMAKE_VERSION VERSION_GREATER 3.1.0 OR CMAKE_VERSION VERSION_EQUAL 3.1.0) # only valid for new platforms set(CMAKE_CXX_STANDARD 11) else() if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang" OR${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") # for older cmake versions # (note, this CXX flag is only valid for clang and gcc, for MSVC it is not needed) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") endif() endif() # enable some warnings add_definitions(-pedantic -Wall) # disable some warnings if("${CMAKE_CXX_COMPILER_ID}" MATCHES "Clang") # (this is only applicable on clang) # ignore some mathgl warnings add_definitions(-Wno-return-type-c-linkage -Wno-keyword-macro -Wno-missing-braces) #set(EXTERNAL_PROJECT_CXX_FLAGS ${EXTERNAL_PROJECT_CXX_FLAGS} -Wno-return-type-c-linkage) endif() if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # disable some warnings on gcc if (CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 6) add_definitions(-Wno-ignored-attributes -Wno-misleading-indentation) #set(EXTERNAL_PROJECT_CXX_FLAGS ${EXTERNAL_PROJECT_CXX_FLAGS} -Wno-ignored-attributes -Wno-misleading-indentation) endif() endif() if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") add_definitions(-Wno-deprecated-declarations -Wno-unknown-pragmas) #set(EXTERNAL_PROJECT_CXX_FLAGS ${EXTERNAL_PROJECT_CXX_FLAGS} -Wno-deprecated-declarations) endif() set(CMAKE_MODULE_PATH${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake-modules/") include(ExternalProject) find_package(Eigen3) if (${EIGEN3_FOUND}) include_directories(${EIGEN3_INCLUDE_DIR}) add_custom_target(Eigen) # dependency dummy else() SET(DOWNLOADING_EIGEN ON) # if not found system wide download message("-- Downloading Eigen3") ExternalProject_Add( Eigen URL http://bitbucket.org/eigen/eigen/get/3.2.7.zip SOURCE_DIR${CMAKE_CURRENT_BINARY_DIR}/Eigen INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/Eigen_install DOWNLOAD_NO_PROGRESS 1 CMAKE_ARGS${EXTERNAL_PROJECT_CMAKE_ARGS_PREFIX} -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/Eigen_install) include_directories(${CMAKE_CURRENT_BINARY_DIR}/Eigen_install/include/eigen3) endif() find_package(MathGL2) if ("${CMAKE_CXX_COMPILER_ID}" MATCHES "Clang") # If mathgl is installed, we do not trust that it is build using clang (in 99% of the cases, # it is not). Therefore, we *always* build mathgl when we use clang set(MATHGL2_FOUND OFF) endif() if (${MATHGL2_FOUND}) add_custom_target(MathGL) # dependency dummy # patch mgl2/config.h file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/mathgl_patched_headers/mgl2) file(READ${MATHGL2_INCLUDE_DIRS}/mgl2/config.h MATHGL_PATCHED_CONFIG_H) STRING(REGEX REPLACE "#define MGL_HAVE_TYPEOF[ \t]+1" "#define MGL_HAVE_TYPEOF 0 // patched" MATHGL_PATCHED_CONFIG_H ${MATHGL_PATCHED_CONFIG_H}) STRING(REGEX REPLACE "#define MGL_HAVE_C99_COMPLEX[ \t]+1" "#define MGL_HAVE_C99_COMPLEX 0 // patched" MATHGL_PATCHED_CONFIG_H${MATHGL_PATCHED_CONFIG_H}) file(WRITE ${PROJECT_BINARY_DIR}/mathgl_patched_headers/mgl2/config.h${MATHGL_PATCHED_CONFIG_H} ) include_directories(${PROJECT_BINARY_DIR}/mathgl_patched_headers/) # use patched config.h for mathgl else() set(DOWNLOADING_MGL ON) message("-- Downloading MathGl") if(NOT EXISTS${CMAKE_CURRENT_BINARY_DIR}/mathgl_install) ExternalProject_Add( MathGL URL http://downloads.sourceforge.net/mathgl/mathgl-2.3.3.tar.gz SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/mathgl_source BINARY_DIR${CMAKE_CURRENT_BINARY_DIR}/mathgl_binary DOWNLOAD_NO_PROGRESS 1 CMAKE_ARGS ${EXTERNAL_PROJECT_CMAKE_ARGS_PREFIX} -DCMAKE_CXX_STANDARD=11 -Denable-openmp=OFF -DMGL_HAVE_TYPEOF=0 -DMGL_HAVE_C99_COMPLEX=0 -DMGL_LIB_INSTALL_DIR=${CMAKE_CURRENT_BINARY_DIR}/mathgl_install/lib/ -DMGL_CGI_PATH=${CMAKE_CURRENT_BINARY_DIR}/mathgl_install/share/mathgl -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/mathgl_install -DCMAKE_INSTALL_NAME_DIR=@rpath INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/mathgl_install ) else() add_custom_target(MathGL) message("Skipped download of mathgl because its already present. Delete ${CMAKE_CURRENT_BINARY_DIR}/mathgl_install to redownload it.") endif() set(MATHGL2_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR}/mathgl_install/include) set(MATHGL2_LIBRARY_NAME${CMAKE_SHARED_LIBRARY_PREFIX}mgl${CMAKE_SHARED_LIBRARY_SUFFIX}) set(MATHGL2_LIBRARIES "${CMAKE_CURRENT_BINARY_DIR}/mathgl_install/lib/${MATHGL2_LIBRARY_NAME}") endif() include_directories(${MATHGL2_INCLUDE_DIRS}) find_package(Figure QUIET) if(FIGURE_FOUND) set(DIRS ${DIRS}${FIGURE_INCLUDE_DIR}) message(STATUS "Function GET_MODULES: Included Figure directory in variable DIRS") set(LIBS ${LIBS}${FIGURE_LIBRARY}) message(STATUS "Function GET_MODULES: Included Figure library in variable LIBS") # case if Figure is not found by FindFigure.cmake - try to get it from MathGL/FigureClass else() set(FIGURE_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/figure) # directory which should contains the source files message(STATUS "Trying to get it from${FIGURE_INCLUDE_DIR} ...") # check if necessary files exist in MathGL/FigureClass foreach(FIGURE_FILE ${FIGURE_FILE_LIST}) if (NOT EXISTS${FIGURE_INCLUDE_DIR}/${FIGURE_FILE}) message(FATAL_ERROR "Could not find necessary files to build Figure library! Try cloning the git repo again or contact someone.") endif() endforeach() message(STATUS "Found necessary Figure files:${FIGURE_INCLUDE_DIR}") add_library(Figure STATIC ${FIGURE_INCLUDE_DIR}/figure.cpp) #add_library(Figure SHARED${FIGURE_INCLUDE_DIR}/figure.cpp) #target_link_libraries(Figure ${MATHGL2_LIBRARIES}) add_dependencies(Figure Eigen) add_dependencies(Figure MathGL) # as libFigure.a was not built yet (this happens when '$ make' is executed) we need to call # target_link_libraries(main Figure) and *not* target_link_libraries(main libFigure.a) set(DIRS ${DIRS}${FIGURE_INCLUDE_DIR}) set(LIBS ${LIBS} Figure) endif() include_directories(${FIGURE_INCLUDE_DIR}) include_directories(${FIGURE_INCLUDE_DIR}/..) ### BlockLU # make blockLU set(SRCS blockLU.cpp) add_executable(blockLU${SRCS}) target_link_libraries(blockLU Figure ${MATHGL2_LIBRARIES})  #include #include #include using namespace Eigen; /* @brief Solve the system Ry=c * for the upper triangular matrix R * This could help you in your implementation * of solve_LSE() * \param[in] R nxn regular, upper triangular matrix * \param[in] c n dim vector * \return y n dim result vector */ /* SAM_LISTING_BEGIN_2 */ VectorXd solve_R(const MatrixXd& R, const VectorXd& c) { int n = R.rows(); assert(n == R.cols() && n == c.size() && "Input dimensions must agree"); // Initialize VectorXd y(n); // Since R is upper triangular, we can solve by backwards substitution for (int i = n-1; i >= 0; --i) { y(i) = c(i); for (int j = n-1; j > i ; --j) { y(i) -= R(i,j) * y(j); } y(i) /= R(i,i); } return y; } /* SAM_LISTING_END_2 */ /* @brief Solve the System Ax=b * for A << R, v, * u.transpose(), 0; * \param[in] R nxn regular, upper triangular matrix * \param[in] v n dim vector * \param[in] u n dim vector * \param[in] b n+1 dim vector * \return x n+1 dim result vector */ /* SAM_LISTING_BEGIN_1 */ VectorXd solve_LSE(const MatrixXd& R, const VectorXd& v, const VectorXd& u, const VectorXd& b) { unsigned n = R.rows(); assert(R.cols() == n && "R has to be square"); assert(n == v.size() && n == u.size() && n+1 == b.size() && "Input dimensions must agree"); // Initialize VectorXd y(n+1), x(n+1); // Solve the system Ax=b by LU-Decomposition. // Solve Ly = b through forward substitution. // Due to the special structure of our L, // the first n entries of y are easy: y.head(n) = b.head(n); // The last element of y is given by$y_n = b_n - u^T\mathbf{R}^{-1}y_{0...n-1}$y(n) = b(n) - u.transpose() * solve_R(R, y.head(n)); // Solve Ux = y by backward substitution. // First we build U MatrixXd U(n+1,n+1); U << R, v, VectorXd::Zero(n).transpose(), -u.transpose()*solve_R(R,v); // Then we solve Ux = y x = solve_R(U,y); // \iffalse Latex comment-delimiter so this part doesn't appear // in the solution and screws up formatting // Note this could be done using less memory by not constructing // U, doing the first step of the back substitution "by hand" and then // calling solve_R() // VectorXd x(n+1); // x(n) = y(n) / (-u.transpose()*solve_R(R,v)); // x.head(n) = solve_R(R,y.head(n) - v*x(n)); //\fi return x; } /* SAM_LISTING_END_1 */ int main() { // Vectors for testing unsigned n = 10; VectorXd v,u,b; u = v = VectorXd::Random(n); b = VectorXd::Random(n+1); // Upper triangular matrix MatrixXd R(n,n); for (unsigned i = 0; i < n; ++i) { for (unsigned j = i; j < n; ++j) { R(i,j) = rand(); //Bad RNG, but sufficient here } } R /= RAND_MAX; //"norm" R for numerical stability // Build matrix A for Eigensolver MatrixXd A(n+1,n+1); A << R, v, u.transpose(), 0; double error = (solve_LSE(R,v,u,b) - A.colPivHouseholderQr().solve(b)).norm(); if (error > 1e-8) { std::cout << "solve_LSE() returns a different solution than Eigen" << std::endl; } else { std::cout << "solve_LSE() and Eigen get the same result" << std::endl; } }  # - Try to find Eigen3 lib # # This module supports requiring a minimum version, e.g. you can do # find_package(Eigen3 3.1.2) # to require version 3.1.2 or newer of Eigen3. # # Once done this will define # # EIGEN3_FOUND - system has eigen lib with correct version # EIGEN3_INCLUDE_DIR - the eigen include directory # EIGEN3_VERSION - eigen version # # This module reads hints about search locations from # the following enviroment variables: # # EIGEN3_ROOT # EIGEN3_ROOT_DIR # Copyright (c) 2006, 2007 Montel Laurent, # Copyright (c) 2008, 2009 Gael Guennebaud, # Copyright (c) 2009 Benoit Jacob # Redistribution and use is allowed according to the terms of the 2-clause BSD license. if(NOT Eigen3_FIND_VERSION) if(NOT Eigen3_FIND_VERSION_MAJOR) set(Eigen3_FIND_VERSION_MAJOR 2) endif(NOT Eigen3_FIND_VERSION_MAJOR) if(NOT Eigen3_FIND_VERSION_MINOR) set(Eigen3_FIND_VERSION_MINOR 91) endif(NOT Eigen3_FIND_VERSION_MINOR) if(NOT Eigen3_FIND_VERSION_PATCH) set(Eigen3_FIND_VERSION_PATCH 0) endif(NOT Eigen3_FIND_VERSION_PATCH) set(Eigen3_FIND_VERSION "${Eigen3_FIND_VERSION_MAJOR}.${Eigen3_FIND_VERSION_MINOR}.${Eigen3_FIND_VERSION_PATCH}") endif(NOT Eigen3_FIND_VERSION) macro(_eigen3_check_version) file(READ "${EIGEN3_INCLUDE_DIR}/Eigen/src/Core/util/Macros.h" _eigen3_version_header) string(REGEX MATCH "define[ \t]+EIGEN_WORLD_VERSION[ \t]+([0-9]+)" _eigen3_world_version_match "${_eigen3_version_header}") set(EIGEN3_WORLD_VERSION "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+EIGEN_MAJOR_VERSION[ \t]+([0-9]+)" _eigen3_major_version_match "${_eigen3_version_header}") set(EIGEN3_MAJOR_VERSION "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+EIGEN_MINOR_VERSION[ \t]+([0-9]+)" _eigen3_minor_version_match "${_eigen3_version_header}") set(EIGEN3_MINOR_VERSION "${CMAKE_MATCH_1}") set(EIGEN3_VERSION${EIGEN3_WORLD_VERSION}.${EIGEN3_MAJOR_VERSION}.${EIGEN3_MINOR_VERSION}) if(${EIGEN3_VERSION} VERSION_LESS${Eigen3_FIND_VERSION}) set(EIGEN3_VERSION_OK FALSE) else(${EIGEN3_VERSION} VERSION_LESS${Eigen3_FIND_VERSION}) set(EIGEN3_VERSION_OK TRUE) endif(${EIGEN3_VERSION} VERSION_LESS${Eigen3_FIND_VERSION}) if(NOT EIGEN3_VERSION_OK) message(STATUS "Eigen3 version ${EIGEN3_VERSION} found in${EIGEN3_INCLUDE_DIR}, " "but at least version ${Eigen3_FIND_VERSION} is required") endif(NOT EIGEN3_VERSION_OK) endmacro(_eigen3_check_version) if (EIGEN3_INCLUDE_DIR) # in cache already _eigen3_check_version() set(EIGEN3_FOUND${EIGEN3_VERSION_OK}) else (EIGEN3_INCLUDE_DIR) # search for signature, 1st try find_path(EIGEN3_INCLUDE_DIR NAMES signature_of_eigen3_matrix_library HINTS ENV EIGEN3_ROOT ENV EIGEN3_ROOT_DIR PATHS ${CMAKE_INSTALL_PREFIX}/include${KDE4_INCLUDE_DIR} PATH_SUFFIXES eigen3 eigen ) if (NOT DEFINED EIGEN3_INCLUDE_DIR)