From 25b7fd9d74e3b22be662828b574c4177017f4ce4 Mon Sep 17 00:00:00 2001 From: zoziha Date: Sat, 6 May 2023 16:37:01 +0800 Subject: [PATCH 1/3] Add solve based on OpenBLAS --- CMakeLists.txt | 13 +++++ README.md | 1 + doc/specs/stdlib_linalg.md | 51 +++++++++++++++++++ example/linalg/CMakeLists.txt | 4 ++ example/linalg/example_solve.f90 | 9 ++++ src/BLAS_common.fypp | 6 +++ src/CMakeLists.txt | 7 +++ src/stdlib_linalg.fypp | 24 +++++++++ src/stdlib_linalg_solve.fypp | 20 ++++++++ test/linalg/CMakeLists.txt | 12 +++++ test/linalg/test_linalg_openblas.fypp | 71 +++++++++++++++++++++++++++ 11 files changed, 218 insertions(+) create mode 100644 example/linalg/example_solve.f90 create mode 100644 src/BLAS_common.fypp create mode 100644 src/stdlib_linalg_solve.fypp create mode 100644 test/linalg/test_linalg_openblas.fypp diff --git a/CMakeLists.txt b/CMakeLists.txt index 471ce6eda..61d123fd7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,6 +67,19 @@ list( "-DPROJECT_VERSION_PATCH=${PROJECT_VERSION_PATCH}" ) +if(NOT DEFINED BUILD_OPENBLAS) + set(BUILD_OPENBLAS OFF CACHE BOOL "Use routines based on OpenBLAS") +endif() +if(BUILD_OPENBLAS) + find_package(PkgConfig REQUIRED) + pkg_check_modules(OPENBLAS REQUIRED openblas) + if(OPENBLAS_FOUND) + list(APPEND fyppFlags "-DUSE_OPENBLAS") + include_directories(${OPENBLAS_INCLUDE_DIRS}) + link_directories(${OPENBLAS_LIBRARY_DIRS}) + endif() +endif() + add_subdirectory(src) if(BUILD_TESTING) diff --git a/README.md b/README.md index 4d912394d..d6df6fb3c 100644 --- a/README.md +++ b/README.md @@ -133,6 +133,7 @@ Important options are Compiling with maximum rank 15 can be resource intensive and requires at least 16 GB of memory to allow parallel compilation or 4 GB memory for sequential compilation. - `-DBUILD_SHARED_LIBS` set to `on` in case you want link your application dynamically against the standard library (default: `off`). - `-DBUILD_TESTING` set to `off` in case you want to disable the stdlib tests (default: `on`). +- `-DBUILD_OPENBLAS` set to `on` in case you want to use routines based on OpenBLAS (default: `off`). - `-DCMAKE_VERBOSE_MAKEFILE` is by default set to `Off`, but if set to `On` will show commands used to compile the code. - `-DCMAKE_BUILD_TYPE` is by default set to `RelWithDebInfo`, which uses compiler flags suitable for code development (but with only `-O2` optimization). Beware the compiler flags set this way will override any compiler flags specified via `FFLAGS`. To prevent this, use `-DCMAKE_BUILD_TYPE=NoConfig` in conjunction with `FFLAGS`. diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 671cfee2f..4fe4b4d6a 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -6,6 +6,17 @@ title: linalg [TOC] +The stdlib linear algebra functions rely on BLAS and LAPACK to provide +efficient low level implementations of standard linear algebra algorithms. + +If possible, highly optimized libraries that take advantage of specialized processor functionality are preferred. +Examples of such libraries are [OpenBLAS][1], [oneAPI MKL(TM)][2]. +Because these libraries are multithreaded and processor dependent, +environment variables may be required to control the number of threads or to specify the processor architecture. + +[1]: https://www.openblas.net/ +[2]: https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html + ## `diag` - Create a diagonal array or extract the diagonal elements of an array ### Status @@ -425,3 +436,43 @@ Specifically, upper Hessenberg matrices satisfy `a_ij = 0` when `j < i-1`, and l ```fortran {!example/linalg/example_is_hessenberg.f90!} ``` + +## `solve` - Solve a-linear-matrix-equation + +### Status + +Experimental + +### Description + +Solve a linear matrix equation, or system of linear scalar equations. + +Note: `solve` supports only single and double precision floating point types. + +### Class + +Impure function. + +### Synatx + +`x = [[stdlib_linalg(module):solve(interface)]](a,b)` + +### Arguments + +`a`: Shall be a rank-2 array of either `real` or `complex` type. This argument is `intent(in)`. Coefficient matrix. + +`b`: Shall be a rank-2 array of either `real` or `complex` type. This argument is `intent(in)`. Ordinate or “dependent variable” values. + +Note: `a` and `b` shall be of the same type. + +### Return value + +Returns a rank-2 array. Solution to the system `ax = b`. Returned shape is identical to `b`. + +### Example + +Solve the system of equations `x1 + 2 * x2 = 1` and `3 * x1 + 5 * x2 = 2`. + +```fortran +{!example/linalg/example_solve.f90!} +``` diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 3f31a5574..9f2cc3652 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -14,3 +14,7 @@ ADD_EXAMPLE(is_symmetric) ADD_EXAMPLE(is_triangular) ADD_EXAMPLE(outer_product) ADD_EXAMPLE(trace) + +if(BUILD_OPENBLAS AND OPENBLAS_FOUND) + ADD_EXAMPLE(solve) +endif() diff --git a/example/linalg/example_solve.f90 b/example/linalg/example_solve.f90 new file mode 100644 index 000000000..43ac30670 --- /dev/null +++ b/example/linalg/example_solve.f90 @@ -0,0 +1,9 @@ +program example_solve + use stdlib_linalg, only: solve + implicit none + real(4) :: a(2, 2), b(2, 1), x(2, 1) + a = reshape([1, 3, 2, 5], [2, 2]) + b = reshape([1, 2], [2, 1]) + x = solve(a, b) + print *, x(:, 1) ! [-1., 1.] +end program example_solve diff --git a/src/BLAS_common.fypp b/src/BLAS_common.fypp new file mode 100644 index 000000000..deef860fc --- /dev/null +++ b/src/BLAS_common.fypp @@ -0,0 +1,6 @@ +#! Define BLAS kinds and types +#:set BLAS_REAL_KINDS = ["sp", "dp"] +#:set BLAS_CMPLX_KINDS = ["sp", "dp"] + +#:set BLAS_REAL_GESV = ["sgesv", "dgesv"] +#:set BLAS_CMPLX_GESV = ["cgesv", "zgesv"] diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ceb1bd2b9..c724bbc7f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,6 +62,12 @@ set(fppFiles stdlib_version.fypp ) +if(BUILD_OPENBLAS AND OPENBLAS_FOUND) + list( + APPEND fppFiles + stdlib_linalg_solve.fypp + ) +endif() fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -85,6 +91,7 @@ set(SRC ) add_library(${PROJECT_NAME} ${SRC}) +target_link_libraries(${PROJECT_NAME} ${OPENBLAS_LIBRARIES}) set_target_properties( ${PROJECT_NAME} diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 3ed905c56..0cf7db17f 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -1,4 +1,5 @@ #:include "common.fypp" +#:include "BLAS_common.fypp" #:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES module stdlib_linalg !!Provides a support for various linear algebra procedures @@ -24,6 +25,10 @@ module stdlib_linalg public :: is_triangular public :: is_hessenberg + #:if defined("USE_OPENBLAS") + public :: solve + #:endif + interface diag !! version: experimental !! @@ -214,6 +219,25 @@ module stdlib_linalg #:endfor end interface is_hessenberg + #:if defined("USE_OPENBLAS") + interface solve + !! version: experimental + !! + !! Solve a linear matrix equation, or system of linear scalar equations + !! ([Specification](../page/specs/stdlib_linalg.html# + !! solve-solve-a-linear-matrix-equation)) + #:set REAL_KINDS_TYPES_GESV = list(zip(BLAS_REAL_KINDS, REAL_TYPES, BLAS_REAL_GESV)) + #:set CMPLX_KINDS_TYPES_GESV = list(zip(BLAS_CMPLX_KINDS, CMPLX_TYPES, BLAS_CMPLX_GESV)) + #:set RC_KINDS_TYPES_GESV = REAL_KINDS_TYPES_GESV + CMPLX_KINDS_TYPES_GESV + #:for k1, t1, g1 in RC_KINDS_TYPES_GESV + module function solve_${t1[0]}$${k1}$(a, b) result(x) + ${t1}$, intent(in) :: a(:, :), b(:, :) + ${t1}$ :: x(size(b, 1), size(b, 2)) + end function solve_${t1[0]}$${k1}$ + #:endfor + end interface solve + #:endif + contains diff --git a/src/stdlib_linalg_solve.fypp b/src/stdlib_linalg_solve.fypp new file mode 100644 index 000000000..05bb3536c --- /dev/null +++ b/src/stdlib_linalg_solve.fypp @@ -0,0 +1,20 @@ +#:include "common.fypp" +#:include "BLAS_common.fypp" +#:set REAL_KINDS_TYPES_GESV = list(zip(BLAS_REAL_KINDS, REAL_TYPES, BLAS_REAL_GESV)) +#:set CMPLX_KINDS_TYPES_GESV = list(zip(BLAS_CMPLX_KINDS, CMPLX_TYPES, BLAS_CMPLX_GESV)) +#:set RC_KINDS_TYPES_GESV = REAL_KINDS_TYPES_GESV + CMPLX_KINDS_TYPES_GESV +submodule (stdlib_linalg) stdlib_linalg_solve + implicit none +contains + #:for k1, t1, g1 in RC_KINDS_TYPES_GESV + module function solve_${t1[0]}$${k1}$(a, b) result(x) + ${t1}$, intent(in) :: a(:, :), b(:, :) + ${t1}$ :: x(size(b, 1), size(b, 2)) + ${t1}$ :: a_local(size(a, 1), size(a, 2)) + integer :: ipiv(size(a, 1)), info + a_local = a + x = b + call ${g1}$(size(a, 1), size(b, 2), a_local, size(a, 1), ipiv, x, size(b, 1), info) + end function solve_${t1[0]}$${k1}$ + #:endfor +end submodule stdlib_linalg_solve diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 4a315f545..a279b5a86 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -3,7 +3,19 @@ set( "test_linalg.fypp" "test_linalg_matrix_property_checks.fypp" ) + +if(BUILD_OPENBLAS AND OPENBLAS_FOUND) + list( + APPEND fppFiles + test_linalg_openblas.fypp + ) +endif() + fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_matrix_property_checks) + +if(BUILD_OPENBLAS AND OPENBLAS_FOUND) + ADDTEST(linalg_openblas) +endif() diff --git a/test/linalg/test_linalg_openblas.fypp b/test/linalg/test_linalg_openblas.fypp new file mode 100644 index 000000000..fd51a9824 --- /dev/null +++ b/test/linalg/test_linalg_openblas.fypp @@ -0,0 +1,71 @@ +#:include "common.fypp" +#:include "BLAS_common.fypp" +#:set REAL_KINDS_TYPES = list(zip(BLAS_REAL_KINDS, REAL_TYPES)) +#:set CMPLX_KINDS_TYPES = list(zip(BLAS_CMPLX_KINDS, CMPLX_TYPES)) +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +module test_linalg_openblas + use testdrive, only: new_unittest, unittest_type, error_type, check, skip_test + use stdlib_kinds, only: sp, dp + use stdlib_linalg, only: solve + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 1000 * epsilon(1._dp) +contains + subroutine collect_linalg_openblas(testsuite) + type(unittest_type), allocatable, intent(out) :: testsuite(:) + #:set IMPLEMENTED_TESTS = ['solve'] + #:set NUM_TESTS = int(len(IMPLEMENTED_TESTS)*len(RC_KINDS_TYPES)) + #! set testsuite dynamically + allocate(testsuite, source=[& + #:set TESTS_WRITTEN = 0 + #:for cur_test in IMPLEMENTED_TESTS + #:for k1, t1 in RC_KINDS_TYPES + #! note that one has to use set directives to increment variable + #:set TESTS_WRITTEN = TESTS_WRITTEN + 1 + #! last test in list should not have comma + #:if TESTS_WRITTEN < NUM_TESTS + new_unittest("${cur_test}$_${t1[0]}$${k1}$", test_${cur_test}$_${t1[0]}$${k1}$), & + #:else + new_unittest("${cur_test}$_${t1[0]}$${k1}$", test_${cur_test}$_${t1[0]}$${k1}$) & + #:endif + #:endfor + #:endfor + ]) + end subroutine collect_linalg_openblas + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_solve_${t1[0]}$${k1}$(error) + type(error_type), allocatable, intent(out) :: error + ${t1}$ :: a(2, 2), b(2, 1), x(2, 1) + a = reshape([1, 3, 2, 5], [2, 2]) + b = reshape([1, 2], [2, 1]) + x = solve(a, b) + call check(error, sum(abs(x(:,1)-[-1, 1])) < ${k1}$tol, "sum(abs(x-[-1, 1])) < ${k1}$tol failed.") + end subroutine test_solve_${t1[0]}$${k1}$ + #:endfor +end module test_linalg_openblas + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_openblas, only : collect_linalg_openblas + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + allocate (testsuites, source=[& + new_testsuite("linalg_openblas", collect_linalg_openblas) & + ]) + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program \ No newline at end of file From a508dc5f8a1eb15b96e0451b204678da1fd01f5f Mon Sep 17 00:00:00 2001 From: zoziha Date: Sat, 6 May 2023 18:09:58 +0800 Subject: [PATCH 2/3] Fix fpm-CI --- ci/fpm-deployment.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/ci/fpm-deployment.sh b/ci/fpm-deployment.sh index bdd3c2b6e..0cac01181 100644 --- a/ci/fpm-deployment.sh +++ b/ci/fpm-deployment.sh @@ -32,6 +32,7 @@ prune=( "$destdir/test/test_hash_functions.f90" "$destdir/src/common.f90" "$destdir/src/f18estop.f90" + "$destdir/example/example_solve.f90" ) major=$(cut -d. -f1 VERSION) From 2aedbcd34786c555c41b819f3502af2c2b3b59c3 Mon Sep 17 00:00:00 2001 From: zoziha Date: Sat, 6 May 2023 20:10:49 +0800 Subject: [PATCH 3/3] Use OpenBLAS in Linux/MINGW CI --- .github/workflows/CI.yml | 33 ++++++++++++++++++++++++++++++ .github/workflows/ci_windows.yml | 35 +++++++++++++++++++++++++++++++- CMakeLists.txt | 4 +--- 3 files changed, 68 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0f53605d2..89ae33b48 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -62,6 +62,12 @@ jobs: brew install gcc@${GCC_V} || brew upgrade gcc@${GCC_V} || true brew link gcc@${GCC_V} + - name: Install OpenBLAS Linux + if: contains( matrix.os, 'ubuntu') + run: | + sudo apt-get install pkg-config + sudo apt-get install libopenblas-dev + - name: Configure with CMake if: ${{ contains(matrix.build, 'cmake') }} run: >- @@ -71,14 +77,32 @@ jobs: -DCMAKE_INSTALL_PREFIX=$PWD/_dist -S . -B ${{ env.BUILD_DIR }} + - name: Configure with CMake and OpenBLAS + if: ${{ contains(matrix.build, 'cmake') && contains( matrix.os, 'ubuntu') }} + run: >- + cmake -Wdev + -DCMAKE_BUILD_TYPE=Release + -DCMAKE_MAXIMUM_RANK:String=4 + -DBUILD_OPENBLAS=ON + -DCMAKE_INSTALL_PREFIX=$PWD/_dist + -S . -B build_with_openblas + - name: Build and compile if: ${{ contains(matrix.build, 'cmake') }} run: cmake --build ${{ env.BUILD_DIR }} --parallel + - name: Build and compile with OpenBLAS + if: ${{ contains(matrix.build, 'cmake') && contains( matrix.os, 'ubuntu') }} + run: cmake --build build_with_openblas --parallel + - name: catch build fail run: cmake --build ${{ env.BUILD_DIR }} --verbose --parallel 1 if: ${{ failure() && contains(matrix.build, 'cmake') }} + - name: catch build fail with OpenBLAS + run: cmake --build build_with_openblas --verbose --parallel 1 + if: ${{ failure() && contains(matrix.build, 'cmake') && contains( matrix.os, 'ubuntu') }} + - name: test if: ${{ contains(matrix.build, 'cmake') }} run: >- @@ -88,6 +112,15 @@ jobs: --output-on-failure --no-tests=error + - name: test with OpenBLAS + if: ${{ contains(matrix.build, 'cmake') && contains( matrix.os, 'ubuntu') }} + run: >- + ctest + --test-dir build_with_openblas + --parallel + --output-on-failure + --no-tests=error + - name: Install project if: ${{ contains(matrix.build, 'cmake') }} run: cmake --install ${{ env.BUILD_DIR }} diff --git a/.github/workflows/ci_windows.yml b/.github/workflows/ci_windows.yml index 1808fea45..8a84c6c1f 100644 --- a/.github/workflows/ci_windows.yml +++ b/.github/workflows/ci_windows.yml @@ -38,6 +38,8 @@ jobs: mingw-w64-${{ matrix.arch }}-python-setuptools mingw-w64-${{ matrix.arch }}-cmake mingw-w64-${{ matrix.arch }}-ninja + mingw-w64-${{ matrix.arch }}-pkg-config + mingw-w64-${{ matrix.arch }}-openblas - name: Setup msys POSIX environment uses: msys2/setup-msys2@v2 @@ -57,7 +59,8 @@ jobs: - name: Install fypp run: pip install fypp - - run: >- + - name: CMake configuration + run: >- cmake -Wdev -B build @@ -88,3 +91,33 @@ jobs: - name: Install project run: cmake --install build + + - name: CMake configuration with OpenBLAS + if: contains(matrix.msystem, 'MINGW') + run: >- + cmake + -Wdev + -B build_with_openblas + -DCMAKE_BUILD_TYPE=Debug + -DCMAKE_Fortran_FLAGS_DEBUG="-Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all -fbacktrace" + -DCMAKE_MAXIMUM_RANK:String=4 + -DBUILD_OPENBLAS=ON + -DCMAKE_INSTALL_PREFIX=$PWD/_dist + + - name: CMake build with OpenBLAS + if: contains(matrix.msystem, 'MINGW') + run: cmake --build build_with_openblas --parallel + + - name: catch build fail with OpenBLAS + if: failure() && contains(matrix.msystem, 'MINGW') + run: cmake --build build_with_openblas --verbose --parallel 1 + + - name: CTest with OpenBLAS + if: contains(matrix.msystem, 'MINGW') + run: ctest --test-dir build_with_openblas --output-on-failure --parallel -V -LE quadruple_precision + + - uses: actions/upload-artifact@v1 + if: failure() && contains(matrix.msystem, 'MINGW') + with: + name: WindowsCMakeTestlog + path: build_with_openblas/Testing/Temporary/LastTest.log diff --git a/CMakeLists.txt b/CMakeLists.txt index 61d123fd7..8d6f1df79 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,9 +67,7 @@ list( "-DPROJECT_VERSION_PATCH=${PROJECT_VERSION_PATCH}" ) -if(NOT DEFINED BUILD_OPENBLAS) - set(BUILD_OPENBLAS OFF CACHE BOOL "Use routines based on OpenBLAS") -endif() +option(BUILD_OPENBLAS "Use routines based on OpenBLAS" OFF) if(BUILD_OPENBLAS) find_package(PkgConfig REQUIRED) pkg_check_modules(OPENBLAS REQUIRED openblas)