Provide BLAS, LAPACK backends and interfaces

.github/workflows/fpm-deployment.yml

@@ -31,9 +31,9 @@ jobs:
         --slave /usr/bin/gcov gcov /usr/bin/gcov-${{ env.GCC_V }}
 
     - name: Install fpm latest release
-      uses: fortran-lang/setup-fpm@v3
+      uses: fortran-lang/setup-fpm@v5
       with:
-        fpm-version: 'v0.4.0'
+        fpm-version: 'v0.10.0'
 
     - name: Run fpm test ⚙
       run: |

ci/fpm.toml

@@ -8,3 +8,6 @@ copyright = "2019-2021 stdlib contributors"
 [dev-dependencies]
 test-drive.git = "https://github.com/fortran-lang/test-drive"
 test-drive.tag = "v0.4.0"
+
+[preprocess]
+[preprocess.cpp]

cmake/stdlib.cmake

@@ -40,3 +40,10 @@ function (fypp_f90 fyppopts fyppfiles f90files)
   preprocess("${FYPP}" "${fyppopts}" "fypp" "f90" "${fyppfiles}" _f90files)
   set(${f90files} ${_f90files} PARENT_SCOPE)
 endfunction()
+
+# For fortran sources that contain C preprocessor flags: create ".F90" files 
+function (fypp_f90pp fyppopts fyppfiles F90files)
+  preprocess("${FYPP}" "${fyppopts}" "fypp" "F90" "${fyppfiles}" _F90files)
+  set(${F90files} ${_F90files} PARENT_SCOPE)
+endfunction()
+

doc/specs/stdlib_linalg.md

@@ -6,6 +6,176 @@ title: linalg
 
 [TOC]
 
+The `stdlib` linear algebra library provides high-level APIs for dealing with common linear algebra operations.
+
+## BLAS and LAPACK
+
+### Status
+
+Experimental
+
+### Description
+
+`BLAS` and `LAPACK` backends provide efficient low level implementations of many linear algebra algorithms, and are employed for non-trivial operators. 
+A Modern Fortran version of the [Reference-LAPACK 3.10.1](http://github.com/reference-LAPACK) implementation is provided as a backend. 
+Modern Fortran modules with full explicit typing features are provided after an 
+[automated conversion](https://github.com/perazz/fortran-lapack/blob/main/scripts/modularize_blas.py) 
+of the legacy codes: 
+- [stdlib_linalg_blas(module)], [stdlib_linalg_lapack(module)] provide kind-agnostic interfaces to all functions.
+- Both libraries are available for 32- (`sp`), 64- (`dp`) and 128-bit (`qp`) `real` and `complex` numbers (the latter if available in the current build)
+- Free format, lower-case style
+- `implicit none(type, external)` applied to all procedures and modules
+- `intent` added and all `pure` procedures where possible
+- All procedure names are prefixed with `stdlib_`, while their generic interface is the same as the BLAS/LAPACK default, with the header character dropped. For example, `stdlib_dgemv`, `stdlib_sgemv`, etc. provide implementations for matrix-vector multiply, while the generic interface is named `gemv`
+- F77-style `parameter`s removed, and all numeric constants have been generalized with KIND-dependent Fortran intrinsics. 
+- preprocessor-based OpenMP directives retained.
+The single-source module structure hopefully allows for cross-procedural inlining which is otherwise impossible without link-time optimization.
+
+When available, highly optimized libraries that take advantage of specialized processor instructions should be preferred over the `stdlib` implementation. 
+Examples of such libraries are: OpenBLAS, MKL (TM), Accelerate, and ATLAS. In order to enable their usage, simply ensure that the following pre-processor macros are defined: 
+
+- `STDLIB_EXTERNAL_BLAS`   wraps all BLAS procedures (except for the 128-bit ones) to an external library
+- `STDLIB_EXTERNAL_LAPACK` wraps all LAPACK procedures (except for the 128-bit ones) to an external library
+
+These can be enabled during the build process. For example, with CMake, one can enable these preprocessor directives using `add_compile_definitions(STDLIB_EXTERNAL_BLAS STDLIB_EXTERNAL_LAPACK)`.
+The same is possible from the `fpm` branch, where the `cpp` preprocessor is enabled by default. For example, the macros can be added to the project's manifest:
+
+```toml
+[dependencies]
+stdlib="*"
+
+# Macros are only needed if using an external library
+[preprocess]
+[preprocess.cpp]
+macros = ["STDLIB_EXTERNAL_BLAS", "STDLIB_EXTERNAL_LAPACK"]
+```
+
+or directly via compiler flags: 
+
+`fpm build --flag "-DSTDLIB_EXTERNAL_BLAS -DSTDLIB_EXTERNAL_LAPACK -lblas -llapack"`.
+
+### Syntax 
+
+All procedures in the `BLAS` and `LAPACK` backends follow the standard interfaces from the 
+[Reference LAPACK](https://www.netlib.org/lapack/). So, the online [Users Guide](https://www.netlib.org/lapack/explore-html/)
+should be consulted for the full API and descriptions of procedure arguments and their usage. 
+
+The `stdlib` implementation makes both kind-agnostic and specific procedure interfaces available via modules
+[stdlib_linalg_blas(module)] and [stdlib_linalg_lapack(module)]. Because all procedures start with a letter 
+[that indicates the base datatype](https://www.netlib.org/lapack/lug/node24.html), the `stdlib` generic
+interface drops the heading letter and contains all kind-dependent implementations. For example, the generic 
+interface to the `axpy` function looks like: 
+
+```fortran  
+!> AXPY: constant times a vector plus a vector.
+interface axpy
+    module procedure stdlib_saxpy
+    module procedure stdlib_daxpy
+    module procedure stdlib_qaxpy
+    module procedure stdlib_caxpy
+    module procedure stdlib_zaxpy
+    module procedure stdlib_waxpy
+end interface axpy
+```
+
+The generic interface is the endpoint for using an external library. Whenever the latter is used, references
+to the internal `module procedure`s are replaced with interfaces to the external library, 
+for example: 
+
+```fortran  
+!> AXPY: constant times a vector plus a vector.
+interface axpy
+    pure subroutine caxpy(n,ca,cx,incx,cy,incy)
+        import sp,dp,qp,ilp,lk 
+        implicit none(type,external) 
+        complex(sp), intent(in) :: ca,cx(*)
+        integer(ilp), intent(in) :: incx,incy,n
+        complex(sp), intent(inout) :: cy(*)
+    end subroutine caxpy
+    ! [....]
+    module procedure stdlib_qaxpy
+end interface axpy
+```
+
+Note that the 128-bit functions are only provided by `stdlib` and always point to the internal implementation. 
+Because 128-bit precision is identified as [stdlib_kinds(module):qp], initials for 128-bit procedures were 
+labelled as `q` (quadruple-precision reals) and `w` ("wide" or quadruple-precision complex numbers). 
+Extended precision ([stdlib_kinds(module):xdp]) calculations are currently not supported.
+
+### Example
+
+```fortran
+{!example/linalg/example_blas_gemv.f90!}
+```
+
+```fortran
+{!example/linalg/example_lapack_getrf.f90!}
+```
+
+### Licensing
+
+The Fortran Standard Library is distributed under the MIT License. `LAPACK` and its contained `BLAS` are a 
+freely-available software package. They are available from [netlib](https://www.netlib.org/lapack/) via anonymous 
+ftp and the World Wide Web. Thus, they can be included in commercial software packages (and have been). 
+The license used for the `BLAS` and `LAPACK` backends is the [modified BSD license](https://www.netlib.org/lapack/LICENSE.txt).
+
+The header of the `LICENSE.txt` file has as its licensing requirements:
+
+    Copyright (c) 1992-2013 The University of Tennessee and The University
+                            of Tennessee Research Foundation.  All rights
+                            reserved.
+    Copyright (c) 2000-2013 The University of California Berkeley. All
+                            rights reserved.
+    Copyright (c) 2006-2013 The University of Colorado Denver.  All rights
+                            reserved.
+
+    $COPYRIGHT$
+
+    Additional copyrights may follow
+
+    $HEADER$
+
+    Redistribution and use in source and binary forms, with or without
+    modification, are permitted provided that the following conditions are
+    met:
+
+    - Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    - Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer listed
+      in this license in the documentation and/or other materials
+      provided with the distribution.
+
+    - Neither the name of the copyright holders nor the names of its
+      contributors may be used to endorse or promote products derived from
+      this software without specific prior written permission.
+
+    The copyright holders provide no reassurances that the source code
+    provided does not infringe any patent, copyright, or any other
+    intellectual property rights of third parties.  The copyright holders
+    disclaim any liability to any recipient for claims brought against
+    recipient by any third party for infringement of that parties
+    intellectual property rights.
+
+    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+So the license for the `LICENSE.txt` code is compatible with the use of
+modified versions of the code in the Fortran Standard Library under the MIT license.    
+Credit for the `BLAS`, `LAPACK` libraries should be given to the [LAPACK authors](https://www.netlib.org/lapack/contributor-list.html).
+According to the original license, we also changed the name of the routines and commented the changes made 
+to the original.
+
 ## `diag` - Create a diagonal array or extract the diagonal elements of an array
 
 ### Status

example/linalg/example_blas_gemv.f90

@@ -0,0 +1,14 @@
+program example_gemv
+  use stdlib_linalg, only: eye
+  use stdlib_linalg_blas, only: sp,gemv
+  implicit none(type,external)
+  real(sp) :: A(2, 2), B(2)
+  B = [1.0,2.0]
+  A = eye(2)
+
+  ! Use legacy BLAS interface 
+  call gemv('No transpose',m=size(A,1),n=size(A,2),alpha=1.0,a=A,lda=size(A,1),x=B,incx=1,beta=0.0,y=B,incy=1)
+
+  print *, B ! returns 1.0 2.0
+
+end program example_gemv

example/linalg/example_lapack_getrf.f90

@@ -0,0 +1,14 @@
+program example_getrf
+  use stdlib_linalg, only: eye
+  use stdlib_linalg_lapack, only: dp,ilp,getrf
+  implicit none(type,external)
+  real(dp) :: A(3, 3)
+  integer(ilp) :: ipiv(3),info
+
+  A = eye(3)
+
+  ! LAPACK matrix factorization interface (overwrite result)
+  call getrf(size(A,1),size(A,2),A,size(A,1),ipiv,info)
+  print *, info ! info==0: Success!
+
+end program example_getrf

include/common.fypp

@@ -87,7 +87,6 @@
 #! Whether Fortran 90 compatible code should be generated
 #:set VERSION90 = defined('VERSION90')
 
-
 #! Ranks to be generated when templates are created
 #:if not defined('MAXRANK')
   #:if VERSION90

src/CMakeLists.txt

@@ -64,8 +64,29 @@ set(fppFiles
     stdlib_version.fypp
 )
 
+# Preprocessed files to contain preprocessor directives -> .F90 
+set(cppFiles
+    stdlib_linalg_constants.fypp
+    stdlib_linalg_blas.fypp
+    stdlib_linalg_blas_aux.fypp
+    stdlib_linalg_blas_s.fypp
+    stdlib_linalg_blas_d.fypp
+    stdlib_linalg_blas_q.fypp
+    stdlib_linalg_blas_c.fypp
+    stdlib_linalg_blas_z.fypp
+    stdlib_linalg_blas_w.fypp
+    stdlib_linalg_lapack.fypp
+    stdlib_linalg_lapack_aux.fypp
+    stdlib_linalg_lapack_s.fypp
+    stdlib_linalg_lapack_d.fypp
+    stdlib_linalg_lapack_q.fypp
+    stdlib_linalg_lapack_c.fypp
+    stdlib_linalg_lapack_z.fypp
+    stdlib_linalg_lapack_w.fypp
+)
 
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
+fypp_f90pp("${fyppFlags}" "${cppFiles}" outPreprocFiles)
 
 set(SRC
     stdlib_ansi.f90
@@ -85,6 +106,7 @@ set(SRC
     stdlib_quadrature_gauss.f90
     stdlib_stringlist_type.f90
     ${outFiles}
+    ${outPreprocFiles}
 )
 
 add_library(${PROJECT_NAME} ${SRC})