Add solve based on OpenBLAS

.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 }}

.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

CMakeLists.txt

@@ -67,6 +67,17 @@ list(
   "-DPROJECT_VERSION_PATCH=${PROJECT_VERSION_PATCH}"
 )
 
+option(BUILD_OPENBLAS "Use routines based on OpenBLAS" OFF)
+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)

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`.
 

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)

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!}
+```

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()

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

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"]

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}

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
 
 

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

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()