Merge pull request #1362 from mgates3/bdsvdx

bdsvdx: fix docs to match code, D and E are not const
Merge pull request #1354 from ACSimon33/coverage_ci
2026-08-21 20:23:06 -06:00 · 2026-08-21 20:53:03 -04:00 · 2026-08-06 08:28:47 -06:00 · 2026-08-06 08:27:40 -06:00 · 2026-08-06 08:21:26 -06:00 · 2026-08-06 15:35:51 +02:00
3848 changed files with 57577 additions and 20291 deletions
@@ -15,22 +15,28 @@ environment:
  global:
    CONDA_INSTALL_LOCN: C:\\Miniconda37-x64
    CTEST_OUTPUT_ON_FAILURE: 1
+  matrix:
+  - BUILD_DEFAULT_API: "ON"
+    BUILD_INDEX64_EXT_API: "OFF"
+  - BUILD_DEFAULT_API: "OFF"
+    BUILD_INDEX64_EXT_API: "ON"

 install:
  - call %CONDA_INSTALL_LOCN%\Scripts\activate.bat
 #  - conda config --set auto_update_conda false
-  - conda install -c conda-forge --yes --quiet flang=11.0.1 jom
-  - call "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" amd64
+  - conda install -c conda-forge --yes --quiet flang flang-rt_win-64 cmake ninja
+  - call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" amd64
  - set "LIB=%CONDA_INSTALL_LOCN%\Library\lib;%LIB%"
  - set "CPATH=%CONDA_INSTALL_LOCN%\Library\include;%CPATH%"

 before_build:
  - ps: if (-Not (Test-Path .\build)) { mkdir build }
  - cd build
-  - cmake -G "NMake Makefiles JOM" -DCMAKE_Fortran_COMPILER=flang -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=ON ..
+  - cmake -G "Ninja" -DCMAKE_Fortran_COMPILER=flang -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=ON -DCBLAS=ON -DLAPACKE=ON -DLAPACKE_WITH_TMG=ON -DBUILD_DEFAULT_API=%BUILD_DEFAULT_API% -DBUILD_INDEX64_EXT_API=%BUILD_INDEX64_EXT_API% ..
+#  - cmake -G "NMake Makefiles JOM" -DCMAKE_Fortran_COMPILER=flang -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=ON ..

 build_script:
  - cmake --build .

 test_script:
-  - ctest -j2
+  - ctest -j2 --output-on-failure 
@@ -0,0 +1,65 @@
+using BinaryBuilder, Pkg
+
+haskey(ENV, "BLAS_LAPACK_RELEASE") || error("The environment variable BLAS_LAPACK_RELEASE is not defined.")
+haskey(ENV, "BLAS_LAPACK_COMMIT") || error("The environment variable BLAS_LAPACK_COMMIT is not defined.")
+haskey(ENV, "BLAS_LAPACK_URL") || error("The environment variable BLAS_LAPACK_URL is not defined.")
+
+name = "blas_lapack"
+version = VersionNumber(ENV["BLAS_LAPACK_RELEASE"])
+
+# Collection of sources required to complete build
+sources = [
+    GitSource(ENV["BLAS_LAPACK_URL"], ENV["BLAS_LAPACK_COMMIT"])
+]
+
+# Bash recipe for building across all platforms
+script = raw"""
+cd ${WORKSPACE}/srcdir/lapack
+
+# FortranCInterface_VERIFY fails on macOS, but it's not actually needed for the current build
+sed -i 's/FortranCInterface_VERIFY/# FortranCInterface_VERIFY/g' ./CBLAS/CMakeLists.txt
+sed -i 's/FortranCInterface_VERIFY/# FortranCInterface_VERIFY/g' ./LAPACKE/include/CMakeLists.txt
+
+mkdir build && cd build
+cmake .. \
+   -DCBLAS=ON \
+   -DLAPACKE=ON \
+   -DCMAKE_INSTALL_PREFIX="$prefix" \
+   -DCMAKE_FIND_ROOT_PATH="$prefix" \
+   -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
+   -DCMAKE_BUILD_TYPE=Release \
+   -DBUILD_SHARED_LIBS=OFF \
+   -DBUILD_INDEX64_EXT_API=OFF \
+   -DTEST_FORTRAN_COMPILER=OFF \
+   -DLAPACKE_WITH_TMG=OFF
+
+make -j${nproc}
+make install
+
+install_license $WORKSPACE/srcdir/lapack/LICENSE
+"""
+
+# These are the platforms we will build for by default, unless further
+# platforms are passed in on the command line
+platforms = supported_platforms()
+platforms = expand_gfortran_versions(platforms)
+
+# The products that we will ensure are always built
+products = [
+    FileProduct("lib/libblas.a", :libblas_a),
+    FileProduct("lib/libcblas.a", :libcblas_a),
+    FileProduct("lib/liblapack.a", :liblapack_a),
+    FileProduct("lib/liblapacke.a", :liblapacke_a),
+    # LibraryProduct("libblas", :libblas),
+    # LibraryProduct("libcblas", :libcblas),
+    # LibraryProduct("liblapack", :liblapack),
+    # LibraryProduct("liblapacke", :liblapacke),
+]
+
+# Dependencies that must be installed before this package can be built
+dependencies = [
+    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
+]
+
+# Build the tarballs, and possibly a `build.jl` as well.
+build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")
@@ -0,0 +1,90 @@
+# Version
+haskey(ENV, "BLAS_LAPACK_RELEASE") || error("The environment variable BLAS_LAPACK_RELEASE is not defined.")
+version = VersionNumber(ENV["BLAS_LAPACK_RELEASE"])
+version2 = ENV["BLAS_LAPACK_RELEASE"]
+package = "blas_lapack"
+
+platforms = [
+   ("aarch64-apple-darwin-libgfortran5"  , "lib", "dylib"),
+#  ("aarch64-linux-gnu-libgfortran3"     , "lib", "so"   ),
+#  ("aarch64-linux-gnu-libgfortran4"     , "lib", "so"   ),
+   ("aarch64-linux-gnu-libgfortran5"     , "lib", "so"   ),
+#  ("aarch64-linux-musl-libgfortran3"    , "lib", "so"   ),
+#  ("aarch64-linux-musl-libgfortran4"    , "lib", "so"   ),
+#  ("aarch64-linux-musl-libgfortran5"    , "lib", "so"   ),
+#  ("powerpc64le-linux-gnu-libgfortran3" , "lib", "so"   ),
+#  ("powerpc64le-linux-gnu-libgfortran4" , "lib", "so"   ),
+#  ("powerpc64le-linux-gnu-libgfortran5" , "lib", "so"   ),
+#  ("x86_64-apple-darwin-libgfortran3"   , "lib", "dylib"),
+#  ("x86_64-apple-darwin-libgfortran4"   , "lib", "dylib"),
+   ("x86_64-apple-darwin-libgfortran5"   , "lib", "dylib"),
+#  ("x86_64-linux-gnu-libgfortran3"      , "lib", "so"   ),
+#  ("x86_64-linux-gnu-libgfortran4"      , "lib", "so"   ),
+   ("x86_64-linux-gnu-libgfortran5"      , "lib", "so"   ),
+#  ("x86_64-linux-musl-libgfortran3"     , "lib", "so"   ),
+#  ("x86_64-linux-musl-libgfortran4"     , "lib", "so"   ),
+#  ("x86_64-linux-musl-libgfortran5"     , "lib", "so"   ),
+#  ("x86_64-unknown-freebsd-libgfortran3", "lib", "so"   ),
+#  ("x86_64-unknown-freebsd-libgfortran4", "lib", "so"   ),
+#  ("x86_64-unknown-freebsd-libgfortran5", "lib", "so"   ),
+#  ("x86_64-w64-mingw32-libgfortran3"    , "bin", "dll"  ),
+#  ("x86_64-w64-mingw32-libgfortran4"    , "bin", "dll"  ),
+   ("x86_64-w64-mingw32-libgfortran5"    , "bin", "dll"  ),
+]
+
+
+for (platform, libdir, ext) in platforms
+
+  tarball_name = "$package.v$version.$platform.tar.gz"
+
+  if isfile("products/$(tarball_name)")
+    # Unzip the tarball generated by BinaryBuilder.jl
+    isdir("products/$platform") && rm("products/$platform", recursive=true)
+    mkdir("products/$platform")
+    run(`tar -xzf products/$(tarball_name) -C products/$platform`)
+
+    if isfile("products/$platform/deps.tar.gz")
+      # Unzip the tarball of the dependencies
+      run(`tar -xzf products/$platform/deps.tar.gz -C products/$platform`)
+
+      # Copy the license of each dependency
+      for folder in readdir("products/$platform/deps/licenses")
+        cp("products/$platform/deps/licenses/$folder", "products/$platform/share/licenses/$folder")
+      end
+      rm("products/$platform/deps/licenses", recursive=true)
+
+      # Copy the shared library of each dependency
+      for file in readdir("products/$platform/deps")
+        cp("products/$platform/deps/$file", "products/$platform/$libdir/$file")
+      end
+
+      # Remove the folder used to unzip the tarball of the dependencies
+      rm("products/$platform/deps", recursive=true)
+      rm("products/$platform/deps.tar.gz", recursive=true)
+    end
+
+    # Create the archives *_binaries
+    isfile("$(package)_binaries.$version2.$platform.tar.gz") && rm("$(package)_binaries.$version2.$platform.tar.gz")
+    isfile("$(package)_binaries.$version2.$platform.zip") && rm("$(package)_binaries.$version2.$platform.zip")
+    cd("products/$platform")
+
+    # Create a folder with the version number of the package
+    mkdir("$(package)_binaries.$version2")
+    for folder in ("include", "share", "lib")
+      cp(folder, "$(package)_binaries.$version2/$folder")
+    end
+
+    cd("$(package)_binaries.$version2")
+    if ext == "dll"
+      run(`zip -r --symlinks ../../../$(package)_binaries.$version2.$platform.zip include share lib`)
+    else
+      run(`tar -czf ../../../$(package)_binaries.$version2.$platform.tar.gz include share lib`)
+    end
+    cd("../../..")
+
+    # Remove the folder used to unzip the tarball generated by BinaryBuilder.jl
+    rm("products/$platform", recursive=true)
+  else
+    @warn("The tarball for the platform $platform was not generated!")
+  end
+end
@@ -7,6 +7,7 @@ on:
    - try-github-actions-for-windows
    paths:
    - .github/workflows/cmake.yml
+    - lapack_testing.py
    - '**CMakeLists.txt'
    - 'BLAS/**'
    - 'CBLAS/**'
@@ -21,6 +22,7 @@ on:
  pull_request:
    paths:
    - .github/workflows/cmake.yml
+    - lapack_testing.py
    - '**CMakeLists.txt'
    - 'BLAS/**'
    - 'CBLAS/**'
@@ -48,7 +50,7 @@ jobs:

  test-install-release:
    # Use GNU compilers
-  
+
    # The CMake configure and build commands are platform agnostic and should work equally
    # well on Windows or Mac.  You can convert this to a matrix build if you need
    # cross-platform coverage.
@@ -62,19 +64,16 @@ jobs:
    strategy:
      fail-fast: true
      matrix:
-        os: [ macos-latest, ubuntu-latest, windows-latest ]
+        os: [ macos-latest, ubuntu-latest, ubuntu-24.04-arm, windows-latest ]
        fflags: [
          "-Wall -Wno-unused-dummy-argument -Wno-unused-variable -Wno-unused-label -Werror=conversion -fimplicit-none -frecursive -fcheck=all",
          "-Wall -Wno-unused-dummy-argument -Wno-unused-variable -Wno-unused-label -Werror=conversion -fimplicit-none -frecursive -fcheck=all -fopenmp" ]
-    
+
    steps:
-    
+
    - name: Checkout LAPACK
      uses: actions/checkout@8e5e7e5ab8b370d6c329ec480221332ada57f0ab # v3.5.2

-    - name: Install ninja-build tool
-      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3
-
    - name: Use GCC-14 on MacOS
      if: ${{ matrix.os == 'macos-latest' }}
      run: >
@@ -90,7 +89,8 @@ jobs:
        -D CMAKE_EXE_LINKER_FLAGS="-Wl,--stack=2097152"

    - name: Configure CMake
-      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using
+      # a single-configuration generator such as make or Ninja.
      # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
      run: >
        cmake -B build -G Ninja
@@ -103,9 +103,9 @@ jobs:
        -D BUILD_SHARED_LIBS:BOOL=ON

    - name: Build
-      # Execute tests defined by the CMake configuration.  
+      # Execute tests defined by the CMake configuration.
      # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
-      run: cmake --build build --config ${{env.BUILD_TYPE}}
+      run: cmake --build build

    - name: Test with OpenMP
      working-directory: ${{github.workspace}}/build
@@ -117,16 +117,53 @@ jobs:
      if: ${{ !contains( matrix.fflags, 'openmp' ) && (matrix.os != 'windows-latest') }}
      run: ctest -C ${{env.BUILD_TYPE}} --schedule-random -j2 --output-on-failure --timeout 100

+    - name: Upload test results
+      id: upload-test-results
+      # Uploaded even when the tests failed; that is when the results
+      # are needed most.  The test summary below links to the artifact.
+      if: ${{ !cancelled() && matrix.os != 'windows-latest' }}
+      uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2
+      with:
+        name: test-results-${{ matrix.os }}-${{ contains( matrix.fflags, 'openmp' ) && 'openmp' || 'no-openmp' }}
+        path: |
+          build/TESTING/testing_results.txt
+          build/lapack_testing_junit.xml
+        if-no-files-found: warn
+        retention-days: 14
+
+    - name: Write test summary
+      if: ${{ !cancelled() && matrix.os != 'windows-latest' }}
+      env:
+        ARTIFACT_URL: ${{ steps.upload-test-results.outputs.artifact-url }}
+      run: |
+        cd build 2>/dev/null || exit 0
+        python3 lapack_testing.py -d TESTING --merge-apis --markdown summary.md || true
+        if [ -f summary.md ]; then
+          cat summary.md >> "$GITHUB_STEP_SUMMARY"
+          if [ -n "$ARTIFACT_URL" ]; then
+            printf '\nThe raw output of every test run (`testing_results.txt`) and a JUnit XML report are in the [test-results artifact](%s).\n' "$ARTIFACT_URL" >> "$GITHUB_STEP_SUMMARY"
+          fi
+        fi
+
    - name: Install
+      # Since we use a single configuration generator, Ninja, there is no need to provide
+      # the '--config ${{env.BUILD_TYPE}}' option for the build step in the 'cmake --build' command.
      run: cmake --build build --target install -j2

-  coverage:
-    runs-on: ubuntu-latest 
+  test-extended-api-only:
+    runs-on: ubuntu-latest
+
    env:
-      BUILD_TYPE: Coverage
-      FFLAGS: "-fopenmp"
-    steps:     
-    
+      BUILD_TYPE: Release
+      FFLAGS: "-Wall -Wno-unused-dummy-argument -Wno-unused-variable -Wno-unused-label -Werror=conversion -fimplicit-none -frecursive -fcheck=all"
+
+    strategy:
+      fail-fast: true
+      matrix:
+        shared_libs: [ OFF, ON ]
+
+    steps:
+
    - name: Checkout LAPACK
      uses: actions/checkout@8e5e7e5ab8b370d6c329ec480221332ada57f0ab # v3.5.2

@@ -134,7 +171,73 @@ jobs:
      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3

    - name: Configure CMake
-      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+      run: >
+        cmake -B build -G Ninja
+        -D CMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
+        -D BUILD_SHARED_LIBS:BOOL=${{matrix.shared_libs}}
+        -D BUILD_DEFAULT_API:BOOL=OFF
+        -D BUILD_INDEX64_EXT_API:BOOL=ON
+        -D BUILD_TESTING:BOOL=ON
+        -D CBLAS:BOOL=ON
+        -D LAPACKE:BOOL=ON
+        -D LAPACKE_WITH_TMG:BOOL=ON
+
+    - name: Build
+      run: cmake --build build --config ${{env.BUILD_TYPE}}
+
+    - name: Test
+      working-directory: ${{github.workspace}}/build
+      run: ctest -C ${{env.BUILD_TYPE}} --schedule-random -j2 --output-on-failure --timeout 100
+
+    - name: Upload test results
+      id: upload-test-results
+      if: ${{ !cancelled() }}
+      uses: actions/upload-artifact@ea165f8d65b6e75b540449e92b4886f43607fa02 # v4.6.2
+      with:
+        name: test-results-extended-api-shared-${{ matrix.shared_libs }}
+        path: |
+          build/TESTING/testing_results.txt
+          build/lapack_testing_junit.xml
+        if-no-files-found: warn
+        retention-days: 14
+
+    - name: Write test summary
+      if: ${{ !cancelled() }}
+      env:
+        ARTIFACT_URL: ${{ steps.upload-test-results.outputs.artifact-url }}
+      run: |
+        cd build 2>/dev/null || exit 0
+        python3 lapack_testing.py -d TESTING --merge-apis --markdown summary.md || true
+        if [ -f summary.md ]; then
+          cat summary.md >> "$GITHUB_STEP_SUMMARY"
+          if [ -n "$ARTIFACT_URL" ]; then
+            printf '\nThe raw output of every test run (`testing_results.txt`) and a JUnit XML report are in the [test-results artifact](%s).\n' "$ARTIFACT_URL" >> "$GITHUB_STEP_SUMMARY"
+          fi
+        fi
+
+  coverage:
+
+    runs-on: ubuntu-latest
+
+    env:
+      BUILD_TYPE: Coverage
+      FFLAGS: "-fopenmp"
+
+    steps:
+
+    - name: Checkout LAPACK
+      uses: actions/checkout@8e5e7e5ab8b370d6c329ec480221332ada57f0ab # v3.5.2
+      with:
+        # Codecov cannot determine which commit a report belongs to from a
+        # depth-1 clone of a pull request merge commit.
+        fetch-depth: 2
+
+    - name: Install ninja-build tool
+      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3
+
+    - name: Configure CMake
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using
+      # a single-configuration generator such as make or Ninja.
      # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
      run: >
        cmake -B build -G Ninja
@@ -147,17 +250,75 @@ jobs:
        -D BUILD_SHARED_LIBS:BOOL=ON

    - name: Install
+      # Since we use a single configuration generator, Ninja, there is no need to provide
+      # the '--config ${{env.BUILD_TYPE}}' option for the build step in the 'cmake --build' command.
      run: cmake --build build --target install -j2

+    - name: Test
+      working-directory: ${{github.workspace}}/build
+      run: ctest -C ${{env.BUILD_TYPE}} --schedule-random -j2 --output-on-failure --timeout 1800
+
    - name: Coverage
+      # Since we use a single configuration generator, Ninja, there is no need to provide
+      # the '--config ${{env.BUILD_TYPE}}' option for the build step in the 'cmake --build' command.
+      if: ${{ !cancelled() }}
+      run: cmake --build build --target coverage
+
+    - name: Summarize coverage
+      # The coverage target discards gcov's output, so an entirely empty report
+      # is indistinguishable from a good one unless we look at the numbers.
+      # Print them, and fail if nothing was measured at all.
+      if: ${{ !cancelled() }}
      run: |
-        echo "Coverage"
-        cmake --build build --target coverage
-        bash <(curl -s https://codecov.io/bash) -X gcov
+        gcda=$(find build -name '*.gcda' | wc -l)
+        reports=$(find build -name '*.gcov' | wc -l)
+        echo "counter files (.gcda): ${gcda}"
+        echo "gcov reports (.gcov): ${reports}"
+        if [ "${gcda}" -eq 0 ] || [ "${reports}" -eq 0 ]; then
+          echo "::error::No coverage data was recorded; the report would be empty."
+          exit 1
+        fi
+        # In a gcov report an executed line is prefixed with its execution
+        # count and an unexecuted one with '#####' or '====='; everything else
+        # ('-') is not executable.
+        set -- $(find build -name '*.gcov' -exec cat {} + | awk '
+          /^ *[0-9]+[*]?:/     { hit++; next }
+          /^ *(#####|=====):/  { miss++ }
+          END { printf "%d %d\n", hit + 0, hit + miss + 0 }')
+        hit=$1
+        total=$2
+        if [ "${hit}" -eq 0 ]; then
+          echo "::error::Coverage report is empty: not a single line was executed."
+          exit 1
+        fi
+        percent=$(awk -v h="${hit}" -v t="${total}" 'BEGIN { printf "%.2f", 100 * h / t }')
+        echo "lines executed: ${hit} of ${total} (${percent}%)"
+        printf '## Coverage\n\n%s%% of lines executed (%s of %s) across %s files.\n' \
+          "${percent}" "${hit}" "${total}" "${reports}" >> "$GITHUB_STEP_SUMMARY"
+
+    - name: Upload coverage report to Codecov
+      if: ${{ !cancelled() }}
+      uses: codecov/codecov-action@fb8b3582c8e4def4969c97caa2f19720cb33a72f # v7.0.0
+      with:
+        token: ${{ secrets.CODECOV_TOKEN }}
+        # The .gcov files already exist, so there is no need for the uploader
+        # to run gcov a second time itself.
+        plugins: noop
+        verbose: true
+        # Pull requests from forks have no access to repository secrets, so
+        # they can only upload if the Codecov organization permits tokenless
+        # uploads.  Do not turn that into a CI failure for the contributor.
+        fail_ci_if_error: ${{ secrets.CODECOV_TOKEN != '' }}

  test-install-cblas-lapacke-without-fortran-compiler:
+
    runs-on: ubuntu-latest
+
+    env:
+      BUILD_TYPE: Release
+
    steps:
+
    - name: Checkout LAPACK
      uses: actions/checkout@8e5e7e5ab8b370d6c329ec480221332ada57f0ab # v3.5.2

@@ -171,9 +332,12 @@ jobs:
        sudo apt purge gfortran

    - name: Configure CMake
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using
+      # a single-configuration generator such as make or Ninja.
+      # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
      run: >
        cmake -B build -G Ninja
-        -D CMAKE_BUILD_TYPE=Release
+        -D CMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
        -D CMAKE_INSTALL_PREFIX=${{github.workspace}}/lapack_install
        -D CBLAS:BOOL=ON
        -D LAPACKE:BOOL=ON
@@ -184,10 +348,14 @@ jobs:
        -D BUILD_SHARED_LIBS:BOOL=ON

    - name: Install
+      # Since we use a single configuration generator, Ninja, there is no need to provide
+      # the '--config ${{env.BUILD_TYPE}}' option for the build step in the 'cmake --build' command.
      run: cmake --build build --target install -j2

  memory-check:
+
    runs-on: ubuntu-latest
+
    env:
      BUILD_TYPE: Debug

@@ -198,13 +366,16 @@ jobs:

    - name: Install ninja-build tool
      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3
-    
+
    - name: Install APT packages
      run: |
        sudo apt update
        sudo apt install -y cmake valgrind gfortran

    - name: Configure CMake
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using
+      # a single-configuration generator such as make or Ninja.
+      # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
      run: >
        cmake -B build -G Ninja
        -D CMAKE_BUILD_TYPE=${{env.BUILD_TYPE}}
@@ -216,12 +387,14 @@ jobs:
        -D LAPACK_TESTING_USE_PYTHON:BOOL=OFF

    - name: Build
-      run: cmake --build build --config ${{env.BUILD_TYPE}}
+      # Since we use a single configuration generator, Ninja, there is no need to provide
+      # the '--config ${{env.BUILD_TYPE}}' option for the build step in the 'cmake --build' command.
+      run: cmake --build build

    - name: Test
      working-directory: ${{github.workspace}}/build
      run: |
-        ctest -C ${{env.BUILD_TYPE}} --schedule-random -j2 -T memcheck > memcheck.out
+        ctest -C ${{env.BUILD_TYPE}} --output-on-failure --schedule-random -j2 -T memcheck > memcheck.out
        cat memcheck.out
        if tail -n 1 memcheck.out | grep -q "Memory checking results:"; then
          exit 0
@@ -0,0 +1,239 @@
+name: Release
+
+on:
+  push:
+    # Sequence of patterns matched against refs/tags
+    tags:
+      - 'v*' # Push events to matching v*, i.e. v1.0, v2023.11.15
+
+jobs:
+  build-linux-x64:
+    name: blas / lapack -- Linux (x86_64) -- Release ${{ github.ref_name }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Install Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "1.7"
+          arch: x64
+
+      - name: Set the environment variables BINARYBUILDER_AUTOMATIC_APPLE, BLAS_LAPACK_RELEASE, BLAS_LAPACK_COMMIT
+        shell: bash
+        run: |
+          echo "BINARYBUILDER_AUTOMATIC_APPLE=true" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_RELEASE=${{ github.ref_name }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_COMMIT=${{ github.sha }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_URL=https://github.com/${{ github.repository }}.git" >> $GITHUB_ENV
+
+      - name: Cross-compilation of blas / lapack -- x86_64-linux-gnu-libgfortran5
+        run: |
+          julia --color=no -e 'using Pkg; Pkg.add("BinaryBuilder")'
+          julia --color=no .github/julia/build_tarballs.jl x86_64-linux-gnu-libgfortran5 --verbose
+
+      - name: Archive artifact
+        run: julia --color=no .github/julia/generate_binaries.jl
+
+      - name: Upload artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: blas_lapack_binaries.${{ github.ref_name }}.x86_64-linux-gnu-libgfortran5.tar.gz
+          path: ./blas_lapack_binaries.${{ github.ref_name }}.x86_64-linux-gnu-libgfortran5.tar.gz
+
+  build-linux-aarch64:
+    name: blas / lapack -- Linux (aarch64) -- Release ${{ github.ref_name }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Install Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "1.7"
+          arch: x64
+
+      - name: Set the environment variables BINARYBUILDER_AUTOMATIC_APPLE, BLAS_LAPACK_RELEASE, BLAS_LAPACK_COMMIT
+        shell: bash
+        run: |
+          echo "BINARYBUILDER_AUTOMATIC_APPLE=true" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_RELEASE=${{ github.ref_name }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_COMMIT=${{ github.sha }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_URL=https://github.com/${{ github.repository }}.git" >> $GITHUB_ENV
+
+      - name: Cross-compilation of blas / lapack -- aarch64-linux-gnu-libgfortran5
+        run: |
+          julia --color=no -e 'using Pkg; Pkg.add("BinaryBuilder")'
+          julia --color=no .github/julia/build_tarballs.jl aarch64-linux-gnu-libgfortran5 --verbose
+
+      - name: Archive artifact
+        run: julia --color=no .github/julia/generate_binaries.jl
+
+      - name: Upload artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: blas_lapack_binaries.${{ github.ref_name }}.aarch64-linux-gnu-libgfortran5.tar.gz
+          path: ./blas_lapack_binaries.${{ github.ref_name }}.aarch64-linux-gnu-libgfortran5.tar.gz
+
+  build-windows-x64:
+    name: blas / lapack -- Windows (x86_64) -- Release ${{ github.ref_name }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Install Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "1.7"
+          arch: x64
+
+      - name: Set the environment variables BINARYBUILDER_AUTOMATIC_APPLE, BLAS_LAPACK_RELEASE, BLAS_LAPACK_COMMIT
+        shell: bash
+        run: |
+          echo "BINARYBUILDER_AUTOMATIC_APPLE=true" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_RELEASE=${{ github.ref_name }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_COMMIT=${{ github.sha }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_URL=https://github.com/${{ github.repository }}.git" >> $GITHUB_ENV
+
+      - name: Cross-compilation of blas / lapack -- x86_64-w64-mingw32-libgfortran5
+        run: |
+          julia --color=no -e 'using Pkg; Pkg.add("BinaryBuilder")'
+          julia --color=no .github/julia/build_tarballs.jl x86_64-w64-mingw32-libgfortran5 --verbose
+      - name: Archive artifact
+        run: julia --color=no .github/julia/generate_binaries.jl
+
+      - name: Upload artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: blas_lapack_binaries.${{ github.ref_name }}.x86_64-w64-mingw32-libgfortran5.zip
+          path: ./blas_lapack_binaries.${{ github.ref_name }}.x86_64-w64-mingw32-libgfortran5.zip
+
+  build-mac-x64:
+    name: blas / lapack -- macOS (x86_64) -- Release ${{ github.ref_name }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Install Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "1.7"
+          arch: x64
+
+      - name: Set the environment variables BINARYBUILDER_AUTOMATIC_APPLE, BLAS_LAPACK_RELEASE, BLAS_LAPACK_COMMIT
+        shell: bash
+        run: |
+          echo "BINARYBUILDER_AUTOMATIC_APPLE=true" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_RELEASE=${{ github.ref_name }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_COMMIT=${{ github.sha }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_URL=https://github.com/${{ github.repository }}.git" >> $GITHUB_ENV
+
+      - name: Cross-compilation of blas / lapack -- x86_64-apple-darwin-libgfortran5
+        run: |
+          julia --color=no -e 'using Pkg; Pkg.add("BinaryBuilder")'
+          julia --color=no .github/julia/build_tarballs.jl x86_64-apple-darwin-libgfortran5 --verbose
+
+      - name: Archive artifact
+        run: julia --color=no .github/julia/generate_binaries.jl
+
+      - name: Upload artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: blas_lapack_binaries.${{ github.ref_name }}.x86_64-apple-darwin-libgfortran5.tar.gz
+          path: ./blas_lapack_binaries.${{ github.ref_name }}.x86_64-apple-darwin-libgfortran5.tar.gz
+
+  build-mac-aarch64:
+    name: blas / lapack -- macOS (aarch64) -- Release ${{ github.ref_name }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Install Julia
+        uses: julia-actions/setup-julia@v2
+        with:
+          version: "1.7"
+          arch: x64
+
+      - name: Set the environment variables BINARYBUILDER_AUTOMATIC_APPLE, BLAS_LAPACK_RELEASE, BLAS_LAPACK_COMMIT
+        shell: bash
+        run: |
+          echo "BINARYBUILDER_AUTOMATIC_APPLE=true" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_RELEASE=${{ github.ref_name }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_COMMIT=${{ github.sha }}" >> $GITHUB_ENV
+          echo "BLAS_LAPACK_URL=https://github.com/${{ github.repository }}.git" >> $GITHUB_ENV
+
+      - name: Cross-compilation of blas / lapack -- aarch64-apple-darwin-libgfortran5
+        run: |
+          julia --color=no -e 'using Pkg; Pkg.add("BinaryBuilder")'
+          julia --color=no .github/julia/build_tarballs.jl aarch64-apple-darwin-libgfortran5 --verbose
+
+      - name: Archive artifact
+        run: julia --color=no .github/julia/generate_binaries.jl
+
+      - name: Upload artifact
+        uses: actions/upload-artifact@v4
+        with:
+          name: blas_lapack_binaries.${{ github.ref_name }}.aarch64-apple-darwin-libgfortran5.tar.gz
+          path: ./blas_lapack_binaries.${{ github.ref_name }}.aarch64-apple-darwin-libgfortran5.tar.gz
+
+  release:
+    name: Create Release and Upload Binaries
+    needs: [build-windows-x64, build-linux-x64, build-linux-aarch64, build-mac-x64, build-mac-aarch64]
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout lapack
+        uses: actions/checkout@v4
+
+      - name: Download artifacts
+        uses: actions/download-artifact@v4
+        with:
+          path: .
+
+      - name: Create GitHub Release
+        run: |
+          gh release create ${{ github.ref_name }} \
+            --title "${{ github.ref_name }}" \
+            --notes "" \
+            --verify-tag
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Upload Linux (x86_64) artifact
+        run: |
+          gh release upload ${{ github.ref_name }} \
+            blas_lapack_binaries.${{ github.ref_name }}.x86_64-linux-gnu-libgfortran5.tar.gz/blas_lapack_binaries.${{ github.ref_name }}.x86_64-linux-gnu-libgfortran5.tar.gz#blas_lapack.${{ github.ref_name }}.linux.x86_64.tar.gz
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Upload Linux (aarch64) artifact
+        run: |
+          gh release upload ${{ github.ref_name }} \
+            blas_lapack_binaries.${{ github.ref_name }}.aarch64-linux-gnu-libgfortran5.tar.gz/blas_lapack_binaries.${{ github.ref_name }}.aarch64-linux-gnu-libgfortran5.tar.gz#blas_lapack.${{ github.ref_name }}.linux.aarch64.tar.gz
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Upload Mac (x86_64) artifact
+        run: |
+          gh release upload ${{ github.ref_name }} \
+            blas_lapack_binaries.${{ github.ref_name }}.x86_64-apple-darwin-libgfortran5.tar.gz/blas_lapack_binaries.${{ github.ref_name }}.x86_64-apple-darwin-libgfortran5.tar.gz#blas_lapack.${{ github.ref_name }}.mac.x86_64.tar.gz
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Upload Mac (aarch64) artifact
+        run: |
+          gh release upload ${{ github.ref_name }} \
+            blas_lapack_binaries.${{ github.ref_name }}.aarch64-apple-darwin-libgfortran5.tar.gz/blas_lapack_binaries.${{ github.ref_name }}.aarch64-apple-darwin-libgfortran5.tar.gz#blas_lapack.${{ github.ref_name }}.mac.aarch64.tar.gz
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Upload Windows (x86_64) artifact
+        run: |
+          gh release upload ${{ github.ref_name }} \
+            blas_lapack_binaries.${{ github.ref_name }}.x86_64-w64-mingw32-libgfortran5.zip/blas_lapack_binaries.${{ github.ref_name }}.x86_64-w64-mingw32-libgfortran5.zip#blas_lapack.${{ github.ref_name }}.windows.x86_64.zip
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
@@ -32,12 +32,12 @@ jobs:

    steps:
      - name: "Checkout code"
-        uses: actions/checkout@c85c95e3d7251135ab7dc9ce3241c5835cc595a9 # v3.5.3
+        uses: actions/checkout@d632683dd7b4114ad314bca15554477dd762a938 # tag=v4.2.0
        with:
          persist-credentials: false

      - name: "Run analysis"
-        uses: ossf/scorecard-action@08b4669551908b1024bb425080c797723083c031 # v2.2.0
+        uses: ossf/scorecard-action@62b2cac7ed8198b15735ed49ab1e5cf35480ba46 # v2.4.0
        with:
          results_file: results.sarif
          results_format: sarif
@@ -59,7 +59,7 @@ jobs:
      # Upload the results as artifacts (optional). Commenting out will disable uploads of run results in SARIF
      # format to the repository Actions tab.
      - name: "Upload artifact"
-        uses: actions/upload-artifact@0b7f8abb1508181956e8e162db84b466c27e18ce # v3.1.2
+        uses: actions/upload-artifact@b4b15b8c7c6ac21ea08fcf65892d2ee8f75cf882 # v4.4.3
        with:
          name: SARIF file
          path: results.sarif
@@ -67,6 +67,6 @@ jobs:

      # Upload the results to GitHub's code scanning dashboard.
      - name: "Upload to code-scanning"
-        uses: github/codeql-action/upload-sarif@f9a7c6738f28efb36e31d49c53a201a9c5d6a476 # v2.14.2
+        uses: github/codeql-action/upload-sarif@662472033e021d55d94146f66f6058822b0b39fd # v3.27.0
        with:
          sarif_file: results.sarif
@@ -1,5 +1,9 @@
 # ignore objects and archives, anywhere in the tree.
 *.[oa]
+*.so
+*.dll
+*.dylib
+*.pdb

 # test in INSTALL
 INSTALL/test*
@@ -23,6 +27,7 @@ CBLAS/examples/cblas_ex1
 CBLAS/examples/cblas_ex2

 # LAPACK testing
+/lapack_testing_junit.xml
 TESTING/LIN/xlintst*
 TESTING/EIG/xeigtst*
 TESTING/EIG/xdmd*
@@ -43,3 +48,7 @@ build*
 DOCS/man
 DOCS/explore-html
 output_err
+
+# Mod files from compilation in SRC
+SRC/la_constants.mod
+SRC/la_xisnan.mod
@@ -29,31 +29,33 @@
 #  Level 1 BLAS
 #---------------------------------------------------------

-set(SBLAS1 isamax.f sasum.f saxpy.f scopy.f sdot.f snrm2.f90
-	srot.f srotg.f90 sscal.f sswap.f sdsdot.f srotmg.f srotm.f)
+set(LAPACK_INSTALL_EXPORT_NAME ${BLASLIB}-targets)

-set(CBLAS1 scabs1.f scasum.f scnrm2.f90 icamax.f caxpy.f ccopy.f
-	cdotc.f cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f)
+set(SBLAS1
+  isamax.f sasum.f saxpy.f saxpby.f scopy.f sdot.f snrm2.f90 srot.f srotg.f90
+  sscal.f sswap.f sdsdot.f srotmg.f srotm.f)

-set(DBLAS1 idamax.f dasum.f daxpy.f dcopy.f ddot.f dnrm2.f90
-	drot.f drotg.f90 dscal.f dsdot.f dswap.f drotmg.f drotm.f)
+set(CBLAS1
+  scabs1.f scasum.f scnrm2.f90 icamax.f90 caxpy.f caxpby.f ccopy.f cdotc.f
+  cdotu.f csscal.f crotg.f90 cscal.f cswap.f csrot.f)
+
+set(DBLAS1
+  idamax.f dasum.f daxpy.f daxpby.f dcopy.f ddot.f dnrm2.f90 drot.f drotg.f90
+  dscal.f dsdot.f dswap.f drotmg.f drotm.f)

 set(DB1AUX sscal.f isamax.f)

-set(ZBLAS1 dcabs1.f dzasum.f dznrm2.f90 izamax.f zaxpy.f zcopy.f
-	zdotc.f zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f)
+set(ZBLAS1
+  dcabs1.f dzasum.f dznrm2.f90 izamax.f90 zaxpy.f zaxpby.f zcopy.f zdotc.f
+  zdotu.f zdscal.f zrotg.f90 zscal.f zswap.f zdrot.f)

 set(CB1AUX
-    isamax.f idamax.f
-    sasum.f saxpy.f scopy.f sdot.f sgemm.f sgemv.f snrm2.f90 srot.f sscal.f
-    sswap.f)
+  isamax.f idamax.f sasum.f saxpy.f scopy.f sdot.f sgemm.f sgemv.f snrm2.f90
+  srot.f sscal.f sswap.f)

 set(ZB1AUX
-    icamax.f idamax.f
-    cgemm.f cherk.f cscal.f ctrsm.f
-    dasum.f daxpy.f dcopy.f ddot.f dgemm.f dgemv.f dnrm2.f90 drot.f dscal.f
-    dswap.f
-    scabs1.f)
+  icamax.f90 idamax.f cgemm.f cherk.f cscal.f ctrsm.f dasum.f daxpy.f dcopy.f
+  ddot.f dgemm.f dgemv.f dnrm2.f90 drot.f dscal.f dswap.f scabs1.f)

 #---------------------------------------------------------------------
 #  Auxiliary routines needed by both the Level 2 and Level 3 BLAS
@@ -63,34 +65,40 @@ set(ALLBLAS lsame.f xerbla.f xerbla_array.f)
 #---------------------------------------------------------
 #  Level 2 BLAS
 #---------------------------------------------------------
-set(SBLAS2 sgemv.f sgbmv.f ssymv.f ssbmv.f sspmv.f
-	strmv.f stbmv.f stpmv.f strsv.f stbsv.f stpsv.f
-	sger.f ssyr.f sspr.f ssyr2.f sspr2.f)
+set(SBLAS2
+  sgemv.f sgbmv.f ssymv.f ssbmv.f sspmv.f strmv.f stbmv.f stpmv.f strsv.f
+  stbsv.f stpsv.f sger.f ssyr.f sspr.f ssyr2.f sspr2.f sskewsymv.f sskewsyr2.f)

-set(CBLAS2 cgemv.f cgbmv.f chemv.f chbmv.f chpmv.f
-	ctrmv.f ctbmv.f ctpmv.f ctrsv.f ctbsv.f ctpsv.f
-	cgerc.f cgeru.f cher.f chpr.f cher2.f chpr2.f)
+set(CBLAS2
+  cgemv.f cgbmv.f chemv.f chbmv.f chpmv.f ctrmv.f ctbmv.f ctpmv.f ctrsv.f
+  ctbsv.f ctpsv.f cgerc.f cgeru.f cher.f chpr.f cher2.f chpr2.f)

-set(DBLAS2 dgemv.f dgbmv.f dsymv.f dsbmv.f dspmv.f
-	dtrmv.f dtbmv.f dtpmv.f dtrsv.f dtbsv.f dtpsv.f
-	dger.f dsyr.f dspr.f dsyr2.f dspr2.f)
+set(DBLAS2
+  dgemv.f dgbmv.f dsymv.f dsbmv.f dspmv.f dtrmv.f dtbmv.f dtpmv.f dtrsv.f
+  dtbsv.f dtpsv.f dger.f dsyr.f dspr.f dsyr2.f dspr2.f dskewsymv.f dskewsyr2.f)

-set(ZBLAS2 zgemv.f zgbmv.f zhemv.f zhbmv.f zhpmv.f
-	ztrmv.f ztbmv.f ztpmv.f ztrsv.f ztbsv.f ztpsv.f
-	zgerc.f zgeru.f zher.f zhpr.f zher2.f zhpr2.f)
+set(ZBLAS2
+  zgemv.f zgbmv.f zhemv.f zhbmv.f zhpmv.f ztrmv.f ztbmv.f ztpmv.f ztrsv.f
+  ztbsv.f ztpsv.f zgerc.f zgeru.f zher.f zhpr.f zher2.f zhpr2.f)

 #---------------------------------------------------------
 #  Level 3 BLAS
 #---------------------------------------------------------
-set(SBLAS3 sgemm.f ssymm.f ssyrk.f ssyr2k.f strmm.f strsm.f sgemmtr.f)
+set(SBLAS3
+  sgemm.f ssymm.f ssyrk.f ssyr2k.f strmm.f strsm.f sgemmtr.f sskewsymm.f
+  sskewsyr2k.f)

-set(CBLAS3 cgemm.f csymm.f csyrk.f csyr2k.f ctrmm.f ctrsm.f
-	chemm.f cherk.f cher2k.f cgemmtr.f)
+set(CBLAS3
+  cgemm.f csymm.f csyrk.f csyr2k.f ctrmm.f ctrsm.f chemm.f cherk.f cher2k.f
+  cgemmtr.f)

-set(DBLAS3 dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.f dgemmtr.f)
+set(DBLAS3
+  dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.f dgemmtr.f dskewsymm.f
+  dskewsyr2k.f)

-set(ZBLAS3 zgemm.f zsymm.f zsyrk.f zsyr2k.f ztrmm.f ztrsm.f
-	zhemm.f zherk.f zher2k.f zgemmtr.f)
+set(ZBLAS3
+  zgemm.f zsymm.f zsyrk.f zsyr2k.f ztrmm.f ztrsm.f zhemm.f zherk.f zher2k.f
+  zgemmtr.f)


 set(SOURCES)
@@ -109,53 +117,43 @@ if(BUILD_COMPLEX16)
 endif()
 list(REMOVE_DUPLICATES SOURCES)

-add_library(${BLASLIB}_obj OBJECT ${SOURCES})
-set_target_properties(${BLASLIB}_obj PROPERTIES POSITION_INDEPENDENT_CODE ON)
+if(BUILD_DEFAULT_API)
+  add_library(${BLASLIB}_obj OBJECT ${SOURCES})
+endif()

 if(BUILD_INDEX64_EXT_API)
-  set(SOURCES_64_F)
-  # Copy files so we can set source property specific to /${BLASLIB}_64_obj target
-  file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${BLASLIB}_64_obj)
-  file(COPY ${SOURCES} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/${BLASLIB}_64_obj)
-  file(GLOB SOURCES_64_F ${CMAKE_CURRENT_BINARY_DIR}/${BLASLIB}_64_obj/*.f*)
-  add_library(${BLASLIB}_64_obj OBJECT ${SOURCES_64_F})
+  include(ExtendedAPIHelpers)
+  generate_64bit_suffixed_sources(${BLASLIB} SOURCES SOURCES_64)
+
+  add_library(${BLASLIB}_64_obj OBJECT ${SOURCES_64})
  target_compile_options(${BLASLIB}_64_obj PRIVATE ${FOPT_ILP64})
-  set_target_properties(${BLASLIB}_64_obj PROPERTIES POSITION_INDEPENDENT_CODE ON)
-  #Add _64 suffix to all Fortran functions via macros
-  foreach(F IN LISTS SOURCES_64_F)
-    if(CMAKE_Fortran_COMPILER_ID STREQUAL "NAG")
-      set_source_files_properties(${F} PROPERTIES COMPILE_FLAGS "-fpp")
-    else()
-      set_source_files_properties(${F} PROPERTIES COMPILE_FLAGS "-cpp")
-    endif()
-    file(STRINGS ${F} ${F}.lst)
-    list(FILTER ${F}.lst INCLUDE REGEX "subroutine|SUBROUTINE|external|EXTERNAL|function|FUNCTION")
-    list(FILTER ${F}.lst EXCLUDE REGEX "^!.*")
-    list(FILTER ${F}.lst EXCLUDE REGEX "^[*].*")
-    list(FILTER ${F}.lst EXCLUDE REGEX "end|END")
-    foreach(FUNC IN LISTS ${F}.lst)
-      string(REGEX REPLACE "^[a-zA-Z0-9_ *]*(subroutine|SUBROUTINE|external|EXTERNAL|function|FUNCTION)[ ]*[*]?" "" FUNC ${FUNC})
-      string(REGEX REPLACE "[(][a-zA-Z0-9_, )]*$" "" FUNC ${FUNC})
-      string(STRIP ${FUNC} FUNC)
-      list(APPEND COPT_64_F "${FUNC}=${FUNC}_64")
-    endforeach()
-    list(REMOVE_DUPLICATES COPT_64_F)
-    set_source_files_properties(${F} PROPERTIES COMPILE_DEFINITIONS "${COPT_64_F}")
-  endforeach()
 endif()

 add_library(${BLASLIB}
-        $<TARGET_OBJECTS:${BLASLIB}_obj>
-        $<$<BOOL:${BUILD_INDEX64_EXT_API}>: $<TARGET_OBJECTS:${BLASLIB}_64_obj>>)
+  $<$<BOOL:${BUILD_DEFAULT_API}>: $<TARGET_OBJECTS:${BLASLIB}_obj>>
+  $<$<BOOL:${BUILD_INDEX64_EXT_API}>: $<TARGET_OBJECTS:${BLASLIB}_64_obj>>)
+
+# For flang, use C linker instead of Fortran linker to avoid macOS-specific flags
+# that CMake adds (tested CMake 4.2).
+if(CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang")
+  set_target_properties (${BLASLIB} PROPERTIES LINKER_LANGUAGE C)
+endif()

 set_target_properties(
  ${BLASLIB} PROPERTIES
  VERSION ${LAPACK_VERSION}
  SOVERSION ${LAPACK_MAJOR_VERSION}
-  POSITION_INDEPENDENT_CODE ON
  )
 lapack_install_library(${BLASLIB})

+add_library(BLAS::BLAS ALIAS ${BLASLIB})
+install(EXPORT ${BLASLIB}-targets
+  FILE ${BLASLIB}-targets.cmake
+  NAMESPACE BLAS::
+  DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/${LAPACKLIB}-${LAPACK_VERSION}
+  COMPONENT Development
+)
+
 if( TEST_FORTRAN_COMPILER )
  add_dependencies( ${BLASLIB} run_test_zcomplexabs run_test_zcomplexdiv run_test_zcomplexmult run_test_zminMax )
 endif()
@@ -69,19 +69,19 @@ all: $(BLASLIB)
 #  Comment out the next 6 definitions if you already have
 #  the Level 1 BLAS.
 #---------------------------------------------------------
-SBLAS1 = isamax.o sasum.o saxpy.o scopy.o sdot.o snrm2.o \
+SBLAS1 = isamax.o sasum.o saxpy.o saxpby.o scopy.o sdot.o snrm2.o \
 	srot.o srotg.o sscal.o sswap.o sdsdot.o srotmg.o srotm.o
 $(SBLAS1): $(FRC)

-CBLAS1 = scabs1.o scasum.o scnrm2.o icamax.o caxpy.o ccopy.o \
+CBLAS1 = scabs1.o scasum.o scnrm2.o icamax.o caxpy.o caxpby.o ccopy.o \
 	cdotc.o cdotu.o csscal.o crotg.o cscal.o cswap.o csrot.o
 $(CBLAS1): $(FRC)

-DBLAS1 = idamax.o dasum.o daxpy.o dcopy.o ddot.o dnrm2.o \
+DBLAS1 = idamax.o dasum.o daxpy.o daxpby.o dcopy.o ddot.o dnrm2.o \
 	drot.o drotg.o dscal.o dsdot.o dswap.o drotmg.o drotm.o
 $(DBLAS1): $(FRC)

-ZBLAS1 = dcabs1.o dzasum.o dznrm2.o izamax.o zaxpy.o zcopy.o \
+ZBLAS1 = dcabs1.o dzasum.o dznrm2.o izamax.o zaxpy.o zaxpby.o zcopy.o \
 	zdotc.o zdotu.o zdscal.o zrotg.o zscal.o zswap.o zdrot.o
 $(ZBLAS1): $(FRC)

@@ -105,7 +105,8 @@ $(ALLBLAS): $(FRC)
 #---------------------------------------------------------
 SBLAS2 = sgemv.o sgbmv.o ssymv.o ssbmv.o sspmv.o \
 	strmv.o stbmv.o stpmv.o strsv.o stbsv.o stpsv.o \
-	sger.o ssyr.o sspr.o ssyr2.o sspr2.o
+	sger.o ssyr.o sspr.o ssyr2.o sspr2.o \
+	sskewsymv.o sskewsyr2.o
 $(SBLAS2): $(FRC)

 CBLAS2 = cgemv.o cgbmv.o chemv.o chbmv.o chpmv.o \
@@ -115,7 +116,8 @@ $(CBLAS2): $(FRC)

 DBLAS2 = dgemv.o dgbmv.o dsymv.o dsbmv.o dspmv.o \
 	dtrmv.o dtbmv.o dtpmv.o dtrsv.o dtbsv.o dtpsv.o \
-	dger.o dsyr.o dspr.o dsyr2.o dspr2.o
+	dger.o dsyr.o dspr.o dsyr2.o dspr2.o \
+	dskewsymv.o dskewsyr2.o
 $(DBLAS2): $(FRC)

 ZBLAS2 = zgemv.o zgbmv.o zhemv.o zhbmv.o zhpmv.o \
@@ -127,14 +129,16 @@ $(ZBLAS2): $(FRC)
 #  Comment out the next 4 definitions if you already have
 #  the Level 3 BLAS.
 #---------------------------------------------------------
-SBLAS3 = sgemm.o ssymm.o ssyrk.o ssyr2k.o strmm.o strsm.o sgemmtr.o
+SBLAS3 = sgemm.o ssymm.o ssyrk.o ssyr2k.o strmm.o strsm.o sgemmtr.o \
+	sskewsymm.o sskewsyr2k.o
 $(SBLAS3): $(FRC)

 CBLAS3 = cgemm.o csymm.o csyrk.o csyr2k.o ctrmm.o ctrsm.o \
 	chemm.o cherk.o cher2k.o cgemmtr.o
 $(CBLAS3): $(FRC)

-DBLAS3 = dgemm.o dsymm.o dsyrk.o dsyr2k.o dtrmm.o dtrsm.o dgemmtr.o
+DBLAS3 = dgemm.o dsymm.o dsyrk.o dsyr2k.o dtrmm.o dtrsm.o dgemmtr.o \
+	dskewsymm.o dskewsyr2k.o
 $(DBLAS3): $(FRC)

 ZBLAS3 = zgemm.o zsymm.o zsyrk.o zsyr2k.o ztrmm.o ztrsm.o \
@@ -0,0 +1,144 @@
+*> \brief \b CAXPBY
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CAXPBY(N,CA,CX,INCX,CB,CY,INCY)
+*
+*       .. Scalar Arguments ..
+*       COMPLEX CA,CB
+*       INTEGER INCX,INCY,N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX CX(*),CY(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*>    CAXPBY constant times a vector plus constant times a vector.
+*>
+*>    Y = ALPHA * X + BETA * Y
+*>
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          number of elements in input vector(s)
+*> \endverbatim
+*>
+*> \param[in] CA
+*> \verbatim
+*>          CA is COMPLEX
+*>          On entry, CA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] CX
+*> \verbatim
+*>          CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          storage spacing between elements of CX
+*> \endverbatim
+*>
+*> \param[in] CB
+*> \verbatim
+*>          CB is COMPLEX
+*>          On entry, CB specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] CY
+*> \verbatim
+*>          CY is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>          storage spacing between elements of CY
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*> \author Martin Koehler, MPI Magdeburg
+*
+*> \ingroup axpby
+*
+*  =====================================================================
+      SUBROUTINE CAXPBY(N,CA,CX,INCX,CB,CY,INCY)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level1 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      COMPLEX CA, CB
+      INTEGER INCX,INCY,N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX CX(*),CY(*)
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL CSCAL
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER I,IX,IY
+*     ..
+      IF (N.LE.0) RETURN
+
+      IF (CA .EQ. (0.0,0.0) .AND. CB.NE.(0.0,0.0)) THEN
+         CALL CSCAL(N,CB, CY, INCY)
+         RETURN
+      END IF
+
+      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
+*
+*        code for both increments equal to 1
+*
+         DO I = 1,N
+            CY(I) = CB*CY(I) + CA*CX(I)
+         END DO
+      ELSE
+*
+*        code for unequal increments or equal increments
+*          not equal to 1
+*
+         IX = 1
+         IY = 1
+         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
+         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
+         DO I = 1,N
+            CY(IY) = CB*CY(IY) + CA*CX(IX)
+            IX = IX + INCX
+            IY = IY + INCY
+         END DO
+      END IF
+*
+      RETURN
+*
+*     End of CAXBPY
+*
+      END
@@ -85,6 +85,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -102,13 +103,16 @@
 *
 *     .. Local Scalars ..
      INTEGER I,IX,IY
+      COMPLEX CDUM
 *     ..
-*     .. External Functions ..
-      REAL SCABS1
-      EXTERNAL SCABS1
+*     .. Statement Functions ..
+      REAL CABS1
+*     ..
+*     .. Statement Function definitions ..
+      CABS1(CDUM) = ABS(REAL(CDUM)) + ABS(AIMAG(CDUM))
 *     ..
      IF (N.LE.0) RETURN
-      IF (SCABS1(CA).EQ.0.0E+0) RETURN
+      IF (CABS1(CA).EQ.0.0E+0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
 *
 *        code for both increments equal to 1
@@ -78,6 +78,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CCOPY(N,CX,INCX,CY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -80,6 +80,7 @@
 *>
 *  =====================================================================
      COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -80,6 +80,7 @@
 *>
 *  =====================================================================
      COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -187,6 +187,7 @@
 *  =====================================================================
      SUBROUTINE CGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,
     +         BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -35,6 +35,16 @@
 *>
 *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
 *> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+*>
+*> Note: if alpha and/or beta is zero, some parts of the matrix-matrix
+*>  operations are not performed. This results in the following NaN/Inf
+*>  propagation quirks:
+*>
+*>  1. If alpha is zero, NaNs or Infs in A or B do not affect the result.
+*>  2. If both alpha and beta are zero, then a zero matrix is returned in C,
+*>   irrespective of any NaNs or Infs in A, B or C.
+*>  3. If only beta is zero, alpha*op( A )*op( B ) is returned, irrespective
+*>   of any NaNs or Infs in C.
 *> \endverbatim
 *
 *  Arguments:
@@ -92,7 +102,9 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is COMPLEX
-*>           On entry, ALPHA specifies the scalar alpha.
+*>           On entry, ALPHA specifies the scalar alpha. If ALPHA is zero the
+*>           values in A and B do not affect the result. This also means that
+*>           NaN/Inf propagation from A and B is inhibited if ALPHA is zero.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -102,7 +114,10 @@
 *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
 *>           part of the array  A  must contain the matrix  A,  otherwise
 *>           the leading  k by m  part of the array  A  must contain  the
-*>           matrix A.
+*>           matrix A, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in A affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of A need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -121,7 +136,10 @@
 *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
 *>           part of the array  B  must contain the matrix  B,  otherwise
 *>           the leading  n by k  part of the array  B  must contain  the
-*>           matrix B.
+*>           matrix B, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in B affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of B need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -136,16 +154,19 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is COMPLEX
-*>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
-*>           supplied as zero then C need not be set on input.
+*>           On entry,  BETA  specifies the scalar  beta. If BETA is zero the
+*>           values in C do not affect the result. This also means that
+*>           NaN/Inf propagation from C is inhibited if BETA is zero.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *> \verbatim
 *>          C is COMPLEX array, dimension ( LDC, N )
 *>           Before entry, the leading  m by n  part of the array  C must
-*>           contain the matrix  C,  except when  beta  is zero, in which
-*>           case C need not be set on entry.
+*>           contain the matrix  C, except if beta is zero.
+*>           If beta is zero, none of the values in C affect the result, even
+*>           if they are NaN/Inf. This also implies that if beta is zero,
+*>           the matrix elements of C need not be initialized by the caller.
 *>           On exit, the array  C  is overwritten by the  m by n  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -185,6 +206,7 @@
 *  =====================================================================
      SUBROUTINE CGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,
     +         BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -50,9 +50,9 @@
 *>           On entry, UPLO specifies whether the lower or the upper
 *>           triangular part of C is access and updated.
 *>
-*>              UPLO = 'L' or 'l', the lower tringular part of C is used.
+*>              UPLO = 'L' or 'l', the lower triangular part of C is used.
 *>
-*>              UPLO = 'U' or 'u', the upper tringular part of C is used.
+*>              UPLO = 'U' or 'u', the upper triangular part of C is used.
 *> \endverbatim
 *
 *> \param[in] TRANSA
@@ -154,7 +154,7 @@
 *>           Before entry, the leading  n by n  part of the array  C must
 *>           contain the matrix  C,  except when  beta  is zero, in which
 *>           case C need not be set on entry.
-*>           On exit, the upper or lower trinangular part of the matrix
+*>           On exit, the upper or lower triangular part of the matrix
 *>           C  is overwritten by the n by n matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -157,6 +157,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -127,6 +127,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CGERC(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -127,6 +127,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CGERU(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -184,6 +184,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -188,6 +188,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -151,6 +151,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHEMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -132,6 +132,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHER(UPLO,N,ALPHA,X,INCX,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -147,6 +147,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHER2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -194,6 +194,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -170,6 +170,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHERK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -146,6 +146,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -127,6 +127,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHPR(UPLO,N,ALPHA,X,INCX,AP)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -142,6 +142,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -86,6 +86,7 @@
 !
 !  =====================================================================
 subroutine CROTG( a, b, c, s )
+   implicit none
   integer, parameter :: wp = kind(1.e0)
 !
 !  -- Reference BLAS level1 routine --
@@ -75,6 +75,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSCAL(N,CA,CX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -95,6 +95,7 @@
 *
 *  =====================================================================
      SUBROUTINE CSROT( N, CX, INCX, CY, INCY, C, S )
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -75,6 +75,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSSCAL(N,SA,CX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -78,6 +78,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSWAP(N,CX,INCX,CY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -186,6 +186,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -185,6 +185,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -164,6 +164,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -183,6 +183,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -197,10 +198,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -271,28 +268,24 @@
              KPLUS1 = K + 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          L = KPLUS1 - J
-                          DO 10 I = MAX(1,J-K),J - 1
-                              X(I) = X(I) + TEMP*A(L+I,J)
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
-                      END IF
+                      TEMP = X(J)
+                      L = KPLUS1 - J
+                      DO 10 I = MAX(1,J-K),J - 1
+                          X(I) = X(I) + TEMP*A(L+I,J)
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          L = KPLUS1 - J
-                          DO 30 I = MAX(1,J-K),J - 1
-                              X(IX) = X(IX) + TEMP*A(L+I,J)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      L = KPLUS1 - J
+                      DO 30 I = MAX(1,J-K),J - 1
+                          X(IX) = X(IX) + TEMP*A(L+I,J)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
                      JX = JX + INCX
                      IF (J.GT.K) KX = KX + INCX
   40             CONTINUE
@@ -300,29 +293,25 @@
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          L = 1 - J
-                          DO 50 I = MIN(N,J+K),J + 1,-1
-                              X(I) = X(I) + TEMP*A(L+I,J)
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(1,J)
-                      END IF
+                      TEMP = X(J)
+                      L = 1 - J
+                      DO 50 I = MIN(N,J+K),J + 1,-1
+                          X(I) = X(I) + TEMP*A(L+I,J)
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(1,J)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          L = 1 - J
-                          DO 70 I = MIN(N,J+K),J + 1,-1
-                              X(IX) = X(IX) + TEMP*A(L+I,J)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(1,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      L = 1 - J
+                      DO 70 I = MIN(N,J+K),J + 1,-1
+                          X(IX) = X(IX) + TEMP*A(L+I,J)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(1,J)
                      JX = JX - INCX
                      IF ((N-J).GE.K) KX = KX - INCX
   80             CONTINUE
@@ -186,6 +186,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -200,10 +201,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -274,59 +271,51 @@
              KPLUS1 = K + 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          L = KPLUS1 - J
-                          IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
-                          TEMP = X(J)
-                          DO 10 I = J - 1,MAX(1,J-K),-1
-                              X(I) = X(I) - TEMP*A(L+I,J)
-   10                     CONTINUE
-                      END IF
+                      L = KPLUS1 - J
+                      IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
+                      TEMP = X(J)
+                      DO 10 I = J - 1,MAX(1,J-K),-1
+                          X(I) = X(I) - TEMP*A(L+I,J)
+   10                 CONTINUE
   20             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 40 J = N,1,-1
                      KX = KX - INCX
-                      IF (X(JX).NE.ZERO) THEN
-                          IX = KX
-                          L = KPLUS1 - J
-                          IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
-                          TEMP = X(JX)
-                          DO 30 I = J - 1,MAX(1,J-K),-1
-                              X(IX) = X(IX) - TEMP*A(L+I,J)
-                              IX = IX - INCX
-   30                     CONTINUE
-                      END IF
+                      IX = KX
+                      L = KPLUS1 - J
+                      IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
+                      TEMP = X(JX)
+                      DO 30 I = J - 1,MAX(1,J-K),-1
+                          X(IX) = X(IX) - TEMP*A(L+I,J)
+                          IX = IX - INCX
+   30                 CONTINUE
                      JX = JX - INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          L = 1 - J
-                          IF (NOUNIT) X(J) = X(J)/A(1,J)
-                          TEMP = X(J)
-                          DO 50 I = J + 1,MIN(N,J+K)
-                              X(I) = X(I) - TEMP*A(L+I,J)
-   50                     CONTINUE
-                      END IF
+                      L = 1 - J
+                      IF (NOUNIT) X(J) = X(J)/A(1,J)
+                      TEMP = X(J)
+                      DO 50 I = J + 1,MIN(N,J+K)
+                          X(I) = X(I) - TEMP*A(L+I,J)
+   50                 CONTINUE
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
                      KX = KX + INCX
-                      IF (X(JX).NE.ZERO) THEN
-                          IX = KX
-                          L = 1 - J
-                          IF (NOUNIT) X(JX) = X(JX)/A(1,J)
-                          TEMP = X(JX)
-                          DO 70 I = J + 1,MIN(N,J+K)
-                              X(IX) = X(IX) - TEMP*A(L+I,J)
-                              IX = IX + INCX
-   70                     CONTINUE
-                      END IF
+                      IX = KX
+                      L = 1 - J
+                      IF (NOUNIT) X(JX) = X(JX)/A(1,J)
+                      TEMP = X(JX)
+                      DO 70 I = J + 1,MIN(N,J+K)
+                          X(IX) = X(IX) - TEMP*A(L+I,J)
+                          IX = IX + INCX
+   70                 CONTINUE
                      JX = JX + INCX
   80             CONTINUE
              END IF
@@ -139,6 +139,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -153,10 +154,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -223,29 +220,25 @@
              KK = 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          K = KK
-                          DO 10 I = 1,J - 1
-                              X(I) = X(I) + TEMP*AP(K)
-                              K = K + 1
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
-                      END IF
+                      TEMP = X(J)
+                      K = KK
+                      DO 10 I = 1,J - 1
+                          X(I) = X(I) + TEMP*AP(K)
+                          K = K + 1
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
                      KK = KK + J
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 30 K = KK,KK + J - 2
-                              X(IX) = X(IX) + TEMP*AP(K)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 30 K = KK,KK + J - 2
+                          X(IX) = X(IX) + TEMP*AP(K)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
                      JX = JX + INCX
                      KK = KK + J
   40             CONTINUE
@@ -254,30 +247,26 @@
              KK = (N* (N+1))/2
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          K = KK
-                          DO 50 I = N,J + 1,-1
-                              X(I) = X(I) + TEMP*AP(K)
-                              K = K - 1
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
-                      END IF
+                      TEMP = X(J)
+                      K = KK
+                      DO 50 I = N,J + 1,-1
+                          X(I) = X(I) + TEMP*AP(K)
+                          K = K - 1
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
                      KK = KK - (N-J+1)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 70 K = KK,KK - (N- (J+1)),-1
-                              X(IX) = X(IX) + TEMP*AP(K)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 70 K = KK,KK - (N- (J+1)),-1
+                          X(IX) = X(IX) + TEMP*AP(K)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
                      JX = JX - INCX
                      KK = KK - (N-J+1)
   80             CONTINUE
@@ -141,6 +141,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -155,10 +156,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -225,29 +222,25 @@
              KK = (N* (N+1))/2
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/AP(KK)
+                      TEMP = X(J)
+                      K = KK - 1
+                      DO 10 I = J - 1,1,-1
+                          X(I) = X(I) - TEMP*AP(K)
+                          K = K - 1
+   10                 CONTINUE
                      KK = KK - J
   20             CONTINUE
              ELSE
                  JX = KX + (N-1)*INCX
                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 30 K = KK - 1,KK - J + 1,-1
+                          IX = IX - INCX
+                          X(IX) = X(IX) - TEMP*AP(K)
+   30                 CONTINUE
                      JX = JX - INCX
                      KK = KK - J
   40             CONTINUE
@@ -256,29 +249,25 @@
              KK = 1
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/AP(KK)
+                      TEMP = X(J)
+                      K = KK + 1
+                      DO 50 I = J + 1,N
+                          X(I) = X(I) - TEMP*AP(K)
+                          K = K + 1
+   50                 CONTINUE
                      KK = KK + (N-J+1)
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 70 K = KK + 1,KK + N - J
+                          IX = IX + INCX
+                          X(IX) = X(IX) - TEMP*AP(K)
+   70                 CONTINUE
                      JX = JX + INCX
                      KK = KK + (N-J+1)
   80             CONTINUE
@@ -174,6 +174,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -275,27 +276,23 @@
              IF (UPPER) THEN
                  DO 50 J = 1,N
                      DO 40 K = 1,M
-                          IF (B(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*B(K,J)
-                              DO 30 I = 1,K - 1
-                                  B(I,J) = B(I,J) + TEMP*A(I,K)
-   30                         CONTINUE
-                              IF (NOUNIT) TEMP = TEMP*A(K,K)
-                              B(K,J) = TEMP
-                          END IF
+                          TEMP = ALPHA*B(K,J)
+                          DO 30 I = 1,K - 1
+                              B(I,J) = B(I,J) + TEMP*A(I,K)
+   30                     CONTINUE
+                          IF (NOUNIT) TEMP = TEMP*A(K,K)
+                          B(K,J) = TEMP
   40                 CONTINUE
   50             CONTINUE
              ELSE
                  DO 80 J = 1,N
                      DO 70 K = M,1,-1
-                          IF (B(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*B(K,J)
-                              B(K,J) = TEMP
-                              IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
-                              DO 60 I = K + 1,M
-                                  B(I,J) = B(I,J) + TEMP*A(I,K)
-   60                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*B(K,J)
+                          B(K,J) = TEMP
+                          IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
+                          DO 60 I = K + 1,M
+                              B(I,J) = B(I,J) + TEMP*A(I,K)
+   60                     CONTINUE
   70                 CONTINUE
   80             CONTINUE
              END IF
@@ -354,12 +351,10 @@
                          B(I,J) = TEMP*B(I,J)
  170                 CONTINUE
                      DO 190 K = 1,J - 1
-                          IF (A(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*A(K,J)
-                              DO 180 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  180                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(K,J)
+                          DO 180 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  180                     CONTINUE
  190                 CONTINUE
  200             CONTINUE
              ELSE
@@ -370,12 +365,10 @@
                          B(I,J) = TEMP*B(I,J)
  210                 CONTINUE
                      DO 230 K = J + 1,N
-                          IF (A(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*A(K,J)
-                              DO 220 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  220                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(K,J)
+                          DO 220 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  220                     CONTINUE
  230                 CONTINUE
  240             CONTINUE
              END IF
@@ -386,16 +379,14 @@
              IF (UPPER) THEN
                  DO 280 K = 1,N
                      DO 260 J = 1,K - 1
-                          IF (A(J,K).NE.ZERO) THEN
-                              IF (NOCONJ) THEN
-                                  TEMP = ALPHA*A(J,K)
-                              ELSE
-                                  TEMP = ALPHA*CONJG(A(J,K))
-                              END IF
-                              DO 250 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  250                         CONTINUE
+                          IF (NOCONJ) THEN
+                              TEMP = ALPHA*A(J,K)
+                          ELSE
+                              TEMP = ALPHA*CONJG(A(J,K))
                          END IF
+                          DO 250 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  250                     CONTINUE
  260                 CONTINUE
                      TEMP = ALPHA
                      IF (NOUNIT) THEN
@@ -414,16 +405,14 @@
              ELSE
                  DO 320 K = N,1,-1
                      DO 300 J = K + 1,N
-                          IF (A(J,K).NE.ZERO) THEN
-                              IF (NOCONJ) THEN
-                                  TEMP = ALPHA*A(J,K)
-                              ELSE
-                                  TEMP = ALPHA*CONJG(A(J,K))
-                              END IF
-                              DO 290 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  290                         CONTINUE
+                          IF (NOCONJ) THEN
+                              TEMP = ALPHA*A(J,K)
+                          ELSE
+                              TEMP = ALPHA*CONJG(A(J,K))
                          END IF
+                          DO 290 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  290                     CONTINUE
  300                 CONTINUE
                      TEMP = ALPHA
                      IF (NOUNIT) THEN
@@ -144,6 +144,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -158,10 +159,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -229,53 +226,45 @@
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          DO 10 I = 1,J - 1
-                              X(I) = X(I) + TEMP*A(I,J)
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(J,J)
-                      END IF
+                      TEMP = X(J)
+                      DO 10 I = 1,J - 1
+                          X(I) = X(I) + TEMP*A(I,J)
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(J,J)
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 30 I = 1,J - 1
-                              X(IX) = X(IX) + TEMP*A(I,J)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 30 I = 1,J - 1
+                          X(IX) = X(IX) + TEMP*A(I,J)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      JX = JX + INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          DO 50 I = N,J + 1,-1
-                              X(I) = X(I) + TEMP*A(I,J)
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(J,J)
-                      END IF
+                      TEMP = X(J)
+                      DO 50 I = N,J + 1,-1
+                          X(I) = X(I) + TEMP*A(I,J)
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(J,J)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 70 I = N,J + 1,-1
-                              X(IX) = X(IX) + TEMP*A(I,J)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 70 I = N,J + 1,-1
+                          X(IX) = X(IX) + TEMP*A(I,J)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      JX = JX - INCX
   80             CONTINUE
              END IF
@@ -177,6 +177,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -209,8 +210,6 @@
      LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER
 *     ..
 *     .. Parameters ..
-      COMPLEX ONE
-      PARAMETER (ONE= (1.0E+0,0.0E+0))
      COMPLEX ZERO
      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
@@ -277,34 +276,26 @@
 *
              IF (UPPER) THEN
                  DO 60 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 30 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-   30                     CONTINUE
-                      END IF
-                      DO 50 K = M,1,-1
-                          IF (B(K,J).NE.ZERO) THEN
-                              IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
-                              DO 40 I = 1,K - 1
-                                  B(I,J) = B(I,J) - B(K,J)*A(I,K)
-   40                         CONTINUE
-                          END IF
+                       DO 30 I = 1,M
+                           B(I,J) = ALPHA*B(I,J)
+   30                  CONTINUE
+                       DO 50 K = M,1,-1
+                           IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+                           DO 40 I = 1,K - 1
+                               B(I,J) = B(I,J) - B(K,J)*A(I,K)
+   40                     CONTINUE
   50                 CONTINUE
   60             CONTINUE
-              ELSE
-                  DO 100 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 70 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-   70                     CONTINUE
-                      END IF
-                      DO 90 K = 1,M
-                          IF (B(K,J).NE.ZERO) THEN
-                              IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
-                              DO 80 I = K + 1,M
-                                  B(I,J) = B(I,J) - B(K,J)*A(I,K)
-   80                         CONTINUE
-                          END IF
+               ELSE
+                   DO 100 J = 1,N
+                       DO 70 I = 1,M
+                           B(I,J) = ALPHA*B(I,J)
+   70                  CONTINUE
+                       DO 90 K = 1,M
+                           IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+                           DO 80 I = K + 1,M
+                               B(I,J) = B(I,J) - B(K,J)*A(I,K)
+   80                     CONTINUE
   90                 CONTINUE
  100             CONTINUE
              END IF
@@ -358,43 +349,33 @@
 *
              IF (UPPER) THEN
                  DO 230 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 190 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-  190                     CONTINUE
-                      END IF
+                      DO 190 I = 1,M
+                          B(I,J) = ALPHA*B(I,J)
+  190                 CONTINUE
                      DO 210 K = 1,J - 1
-                          IF (A(K,J).NE.ZERO) THEN
-                              DO 200 I = 1,M
-                                  B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  200                         CONTINUE
-                          END IF
+                          DO 200 I = 1,M
+                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
+  200                     CONTINUE
  210                 CONTINUE
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(J,J)
                          DO 220 I = 1,M
-                              B(I,J) = TEMP*B(I,J)
+                              B(I,J) = B(I,J)/A(J,J)
  220                     CONTINUE
                      END IF
  230             CONTINUE
              ELSE
                  DO 280 J = N,1,-1
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 240 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-  240                     CONTINUE
-                      END IF
+                      DO 240 I = 1,M
+                          B(I,J) = ALPHA*B(I,J)
+  240                 CONTINUE
                      DO 260 K = J + 1,N
-                          IF (A(K,J).NE.ZERO) THEN
-                              DO 250 I = 1,M
-                                  B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  250                         CONTINUE
-                          END IF
+                          DO 250 I = 1,M
+                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
+  250                     CONTINUE
  260                 CONTINUE
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(J,J)
                          DO 270 I = 1,M
-                              B(I,J) = TEMP*B(I,J)
+                              B(I,J) = B(I,J)/A(J,J)
  270                     CONTINUE
                      END IF
  280             CONTINUE
@@ -408,61 +389,55 @@
                  DO 330 K = N,1,-1
                      IF (NOUNIT) THEN
                          IF (NOCONJ) THEN
-                              TEMP = ONE/A(K,K)
+                              DO 290 I = 1,M
+                                  B(I,K) = B(I,K)/A(K,K)
+  290                         CONTINUE
                          ELSE
-                              TEMP = ONE/CONJG(A(K,K))
+                              DO 390 I = 1,M
+                                  B(I,K) = B(I,K)/CONJG(A(K,K))
+  390                         CONTINUE
                          END IF
-                          DO 290 I = 1,M
-                              B(I,K) = TEMP*B(I,K)
-  290                     CONTINUE
                      END IF
                      DO 310 J = 1,K - 1
-                          IF (A(J,K).NE.ZERO) THEN
-                              IF (NOCONJ) THEN
-                                  TEMP = A(J,K)
-                              ELSE
-                                  TEMP = CONJG(A(J,K))
-                              END IF
-                              DO 300 I = 1,M
-                                  B(I,J) = B(I,J) - TEMP*B(I,K)
-  300                         CONTINUE
+                          IF (NOCONJ) THEN
+                              TEMP = A(J,K)
+                          ELSE
+                              TEMP = CONJG(A(J,K))
                          END IF
+                          DO 300 I = 1,M
+                              B(I,J) = B(I,J) - TEMP*B(I,K)
+  300                     CONTINUE
  310                 CONTINUE
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 320 I = 1,M
-                              B(I,K) = ALPHA*B(I,K)
-  320                     CONTINUE
-                      END IF
+                      DO 320 I = 1,M
+                          B(I,K) = ALPHA*B(I,K)
+  320                 CONTINUE
  330             CONTINUE
              ELSE
                  DO 380 K = 1,N
                      IF (NOUNIT) THEN
                          IF (NOCONJ) THEN
-                              TEMP = ONE/A(K,K)
+                              DO 340 I = 1,M
+                                  B(I,K) = B(I,K)/A(K,K)
+  340                         CONTINUE
                          ELSE
-                              TEMP = ONE/CONJG(A(K,K))
+                              DO 400 I = 1,M
+                                  B(I,K) = B(I,K)/CONJG(A(K,K))
+  400                         CONTINUE
                          END IF
-                          DO 340 I = 1,M
-                              B(I,K) = TEMP*B(I,K)
-  340                     CONTINUE
                      END IF
                      DO 360 J = K + 1,N
-                          IF (A(J,K).NE.ZERO) THEN
-                              IF (NOCONJ) THEN
-                                  TEMP = A(J,K)
-                              ELSE
-                                  TEMP = CONJG(A(J,K))
-                              END IF
-                              DO 350 I = 1,M
-                                  B(I,J) = B(I,J) - TEMP*B(I,K)
-  350                         CONTINUE
+                          IF (NOCONJ) THEN
+                              TEMP = A(J,K)
+                          ELSE
+                              TEMP = CONJG(A(J,K))
                          END IF
+                          DO 350 I = 1,M
+                              B(I,J) = B(I,J) - TEMP*B(I,K)
+  350                     CONTINUE
  360                 CONTINUE
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 370 I = 1,M
-                              B(I,K) = ALPHA*B(I,K)
-  370                     CONTINUE
-                      END IF
+                      DO 370 I = 1,M
+                          B(I,K) = ALPHA*B(I,K)
+  370                 CONTINUE
  380             CONTINUE
              END IF
          END IF
@@ -146,6 +146,7 @@
 *>
 *  =====================================================================
      SUBROUTINE CTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -160,10 +161,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
      COMPLEX TEMP
@@ -231,52 +228,44 @@
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/A(J,J)
-                          TEMP = X(J)
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*A(I,J)
-   10                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/A(J,J)
+                      TEMP = X(J)
+                      DO 10 I = J - 1,1,-1
+                          X(I) = X(I) - TEMP*A(I,J)
+   10                 CONTINUE
   20             CONTINUE
              ELSE
                  JX = KX + (N-1)*INCX
                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/A(J,J)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 I = J - 1,1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*A(I,J)
-   30                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 30 I = J - 1,1,-1
+                          IX = IX - INCX
+                          X(IX) = X(IX) - TEMP*A(I,J)
+   30                 CONTINUE
                      JX = JX - INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/A(J,J)
-                          TEMP = X(J)
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*A(I,J)
-   50                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/A(J,J)
+                      TEMP = X(J)
+                      DO 50 I = J + 1,N
+                          X(I) = X(I) - TEMP*A(I,J)
+   50                 CONTINUE
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/A(J,J)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 I = J + 1,N
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*A(I,J)
-   70                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 70 I = J + 1,N
+                          IX = IX + INCX
+                          X(IX) = X(IX) - TEMP*A(I,J)
+   70                 CONTINUE
                      JX = JX + INCX
   80             CONTINUE
              END IF
@@ -68,6 +68,7 @@
 *>
 *  =====================================================================
      DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -0,0 +1,149 @@
+*> \brief \b DAXPBY
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DAXPBY(N,DA,DX,INCX,DB,DY,INCY)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION DA,DB
+*       INTEGER INCX,INCY,N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION DX(*),DY(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*>    DAXPBY constant times a vector plus constant times a vector.
+*>
+*>    Y = ALPHA * X + BETA * Y
+*>
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          number of elements in input vector(s)
+*> \endverbatim
+*>
+*> \param[in] DA
+*> \verbatim
+*>           DA is DOUBLE PRECISION
+*>           On entry, DA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] DX
+*> \verbatim
+*>          DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          storage spacing between elements of DX
+*> \endverbatim
+*>
+*> \param[in] DB
+*> \verbatim
+*>           DB is DOUBLE PRECISION
+*>           On entry, DB specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] DY
+*> \verbatim
+*>          DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>          storage spacing between elements of DY
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*> \author Martin Koehler, MPI Magdeburg
+*
+*> \ingroup axpby
+*
+*  =====================================================================
+      SUBROUTINE DAXPBY(N,DA,DX,INCX,DB,DY,INCY)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level1 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION DA,DB
+      INTEGER INCX,INCY,N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION DX(*),DY(*)
+*     ..
+*     .. External Subroutines
+      EXTERNAL DSCAL
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER I,IX,IY,M,MP1
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MOD
+*     ..
+      IF (N.LE.0) RETURN
+
+*     Scale if DA.EQ.0
+      IF (DA.EQ.0.0D0 .AND. DB.NE.0.0D0) THEN
+          CALL DSCAL(N, DB, DY, INCY)
+          RETURN
+      END IF
+
+      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
+*
+*        code for both increments equal to 1
+*
+*
+*
+         DO I = 1,N
+            DY(I) = DB*DY(I) + DA*DX(I)
+         END DO
+      ELSE
+*
+*        code for unequal increments or equal increments
+*          not equal to 1
+*
+         IX = 1
+         IY = 1
+         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
+         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
+         DO I = 1,N
+          DY(IY) = DB*DY(IY) + DA*DX(IX)
+          IX = IX + INCX
+          IY = IY + INCY
+         END DO
+      END IF
+      RETURN
+*
+*     End of DAXPBY
+*
+      END
@@ -86,6 +86,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -44,6 +44,7 @@
 *
 *  =====================================================================
      DOUBLE PRECISION FUNCTION DCABS1(Z)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -79,6 +79,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -79,6 +79,7 @@
 *>
 *  =====================================================================
      DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -185,6 +185,7 @@
 *  =====================================================================
      SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,
     +         BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -35,6 +35,16 @@
 *>
 *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
 *> an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
+*>
+*> Note: if alpha and/or beta is zero, some parts of the matrix-matrix
+*>  operations are not performed. This results in the following NaN/Inf
+*>  propagation quirks:
+*>
+*>  1. If alpha is zero, NaNs or Infs in A or B do not affect the result.
+*>  2. If both alpha and beta are zero, then a zero matrix is returned in C,
+*>   irrespective of any NaNs or Infs in A, B or C.
+*>  3. If only beta is zero, alpha*op( A )*op( B ) is returned, irrespective
+*>   of any NaNs or Infs in C.
 *> \endverbatim
 *
 *  Arguments:
@@ -51,6 +61,9 @@
 *>              TRANSA = 'T' or 't',  op( A ) = A**T.
 *>
 *>              TRANSA = 'C' or 'c',  op( A ) = A**T.
+*>
+*>           Note: TRANSA = 'C' is supported for the sake of API consistency
+*>           between all ?GEMM variants.
 *> \endverbatim
 *>
 *> \param[in] TRANSB
@@ -64,6 +77,9 @@
 *>              TRANSB = 'T' or 't',  op( B ) = B**T.
 *>
 *>              TRANSB = 'C' or 'c',  op( B ) = B**T.
+*>
+*>           Note: TRANSB = 'C' is supported for the sake of API consistency
+*>           between all ?GEMM variants.
 *> \endverbatim
 *>
 *> \param[in] M
@@ -92,7 +108,9 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is DOUBLE PRECISION.
-*>           On entry, ALPHA specifies the scalar alpha.
+*>           On entry, ALPHA specifies the scalar alpha. If ALPHA is zero the
+*>           values in A and B do not affect the result. This also means that
+*>           NaN/Inf propagation from A and B is inhibited if ALPHA is zero.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -102,7 +120,10 @@
 *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
 *>           part of the array  A  must contain the matrix  A,  otherwise
 *>           the leading  k by m  part of the array  A  must contain  the
-*>           matrix A.
+*>           matrix A, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in A affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of A need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -121,7 +142,10 @@
 *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
 *>           part of the array  B  must contain the matrix  B,  otherwise
 *>           the leading  n by k  part of the array  B  must contain  the
-*>           matrix B.
+*>           matrix B, except if ALPHA is zero.
+*>           If ALPHA is zero, none of the values in B affect the result, even
+*>           if they are NaN/Inf. This also implies that if ALPHA is zero,
+*>           the matrix elements of B need not be initialized by the caller.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -136,16 +160,19 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is DOUBLE PRECISION.
-*>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
-*>           supplied as zero then C need not be set on input.
+*>           On entry,  BETA  specifies the scalar  beta.  If BETA is zero the
+*>           values in C do not affect the result. This also means that
+*>           NaN/Inf propagation from C is inhibited if BETA is zero.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *> \verbatim
 *>          C is DOUBLE PRECISION array, dimension ( LDC, N )
 *>           Before entry, the leading  m by n  part of the array  C must
-*>           contain the matrix  C,  except when  beta  is zero, in which
-*>           case C need not be set on entry.
+*>           contain the matrix  C, except if beta is zero.
+*>           If beta is zero, none of the values in C affect the result, even
+*>           if they are NaN/Inf. This also implies that if beta is zero,
+*>           the matrix elements of C need not be initialized by the caller.
 *>           On exit, the array  C  is overwritten by the  m by n  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -185,6 +212,7 @@
 *  =====================================================================
      SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,
     +         BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -50,9 +50,9 @@
 *>           On entry, UPLO specifies whether the lower or the upper
 *>           triangular part of C is access and updated.
 *>
-*>              UPLO = 'L' or 'l', the lower tringular part of C is used.
+*>              UPLO = 'L' or 'l', the lower triangular part of C is used.
 *>
-*>              UPLO = 'U' or 'u', the upper tringular part of C is used.
+*>              UPLO = 'U' or 'u', the upper triangular part of C is used.
 *> \endverbatim
 *
 *> \param[in] TRANSA
@@ -154,7 +154,7 @@
 *>           Before entry, the leading  n by n  part of the array  C must
 *>           contain the matrix  C,  except when  beta  is zero, in which
 *>           case C need not be set on entry.
-*>           On exit, the upper or lower trinangular part of the matrix
+*>           On exit, the upper or lower triangular part of the matrix
 *>           C  is overwritten by the n by n matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -426,6 +426,6 @@
 *
      RETURN
 *
-*     End of SGEMM
+*     End of DGEMMTR
 *
      END
@@ -155,6 +155,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -127,6 +127,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -85,11 +85,12 @@
 !> \endverbatim
 !>
 !  =====================================================================
-function DNRM2( n, x, incx ) 
+function DNRM2( n, x, incx )
+   implicit none
   integer, parameter :: wp = kind(1.d0)
   real(wp) :: DNRM2
 !
-!  -- Reference BLAS level1 routine (version 3.9.1) --
+!  -- Reference BLAS level1 routine --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021
@@ -89,6 +89,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DROT(N,DX,INCX,DY,INCY,C,S)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -89,6 +89,7 @@
 !
 !  =====================================================================
 subroutine DROTG( a, b, c, s )
+   implicit none
   integer, parameter :: wp = kind(1.d0)
 !
 !  -- Reference BLAS level1 routine --
@@ -38,6 +38,10 @@
 *>    H=(          )    (          )    (          )    (          )
 *>      (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0).
 *>    SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM.
+*>
+*>    IF DFLAG IS NOT ONE OF THE LISTED ABOVE, THE BEHAVIOR IS UNDEFINED.
+*>    NANS IN DFLAG MAY NOT PROPAGATE TO THE OUTPUT.
+*>
 *> \endverbatim
 *
 *  Arguments:
@@ -93,6 +97,7 @@
 *
 *  =====================================================================
      SUBROUTINE DROTM(N,DX,INCX,DY,INCY,DPARAM)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -24,14 +24,18 @@
 *> \verbatim
 *>
 *>    CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
-*>    THE SECOND COMPONENT OF THE 2-VECTOR  (DSQRT(DD1)*DX1,DSQRT(DD2)*>    DY2)**T.
-*>    WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
+*>    THE SECOND COMPONENT OF THE 2-VECTOR  
+*>         (DSQRT(DD1)*DX1,DSQRT(DD2)*DY2)**T
+*>    WITH DPARAM(1)=DFLAG.
 *>
-*>    DFLAG=-1.D0     DFLAG=0.D0        DFLAG=1.D0     DFLAG=-2.D0
+*>    H HAS ONE OF THE FOLLOWING FORMS:
+*>
+*>      DFLAG=-1.D0     DFLAG=0.D0      DFLAG=1.D0      DFLAG=-2.D0
 *>
 *>      (DH11  DH12)    (1.D0  DH12)    (DH11  1.D0)    (1.D0  0.D0)
 *>    H=(          )    (          )    (          )    (          )
 *>      (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0).
+*>
 *>    LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
 *>    RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
 *>    VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
@@ -87,6 +91,7 @@
 *
 *  =====================================================================
      SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -195,7 +200,7 @@
                  DH11 = ONE
                  DH22 = ONE
                  DFLAG = -ONE
-               ELSE
+               ELSE IF (DFLAG.EQ.ONE) THEN
                  DH21 = -ONE
                  DH12 = ONE
                  DFLAG = -ONE
@@ -220,7 +225,7 @@
                  DH11 = ONE
                  DH22 = ONE
                  DFLAG = -ONE
-               ELSE
+               ELSE IF (DFLAG.EQ.ONE) THEN
                  DH21 = -ONE
                  DH12 = ONE
                  DFLAG = -ONE
@@ -181,6 +181,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -76,6 +76,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSCAL(N,DA,DX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -116,6 +116,7 @@
 *>
 *  =====================================================================
      DOUBLE PRECISION FUNCTION DSDOT(N,SX,INCX,SY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -0,0 +1,365 @@
+*> \brief \b DSKEWSYMM
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DSKEWSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION ALPHA,BETA
+*       INTEGER LDA,LDB,LDC,M,N
+*       CHARACTER SIDE,UPLO
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DSKEWSYMM performs one of the matrix-matrix operations
+*>
+*>    C := alpha*A*B + beta*C,
+*>
+*> or
+*>
+*>    C := alpha*B*A + beta*C,
+*>
+*> where alpha and beta are scalars,  A is a skew-symmetric matrix and  B and
+*> C are  m by n matrices.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>           On entry,  SIDE  specifies whether  the  skew-symmetric matrix  A
+*>           appears on the  left or right  in the  operation as follows:
+*>
+*>              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
+*>
+*>              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
+*> \endverbatim
+*>
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
+*>           triangular  part  of  the  skew-symmetric  matrix   A  is  to  be
+*>           referenced as follows:
+*>
+*>              UPLO = 'U' or 'u'   Only the upper triangular part of the
+*>                                  skew-symmetric matrix is to be referenced.
+*>
+*>              UPLO = 'L' or 'l'   Only the lower triangular part of the
+*>                                  skew-symmetric matrix is to be referenced.
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>           On entry,  M  specifies the number of rows of the matrix  C.
+*>           M  must be at least zero.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>           On entry, N specifies the number of columns of the matrix C.
+*>           N  must be at least zero.
+*> \endverbatim
+*>
+*> \param[in] ALPHA
+*> \verbatim
+*>          ALPHA is DOUBLE PRECISION
+*>           On entry, ALPHA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
+*>           m  when  SIDE = 'L' or 'l'  and is  n otherwise.
+*>           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
+*>           the array  A  must contain the  skew-symmetric matrix,  such that
+*>           when  UPLO = 'U' or 'u', the strictly m by m upper triangular
+*>           part of the array  A  must contain the upper triangular part
+*>           of the  skew-symmetric matrix and the  leading  lower triangular
+*>           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
+*>           the strictly  m by m  lower triangular part  of the  array  A
+*>           must  contain  the  lower triangular part  of the  skew-symmetric
+*>           matrix and the leading  upper triangular part of  A  is not
+*>           referenced.
+*>           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
+*>           the array  A  must contain the  skew-symmetric matrix,  such that
+*>           when  UPLO = 'U' or 'u', the strictly n by n upper triangular
+*>           part of the array  A  must contain the upper triangular part
+*>           of the  skew-symmetric matrix and the  leading  lower triangular
+*>           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
+*>           the strictly  n by n  lower triangular part  of the  array  A
+*>           must  contain  the  lower triangular part  of the  skew-symmetric
+*>           matrix and the leading  upper triangular part of  A  is not
+*>           referenced.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>           On entry, LDA specifies the first dimension of A as declared
+*>           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
+*>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
+*>           least  max( 1, n ).
+*> \endverbatim
+*>
+*> \param[in] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension ( LDB, N )
+*>           Before entry, the leading  m by n part of the array  B  must
+*>           contain the matrix B.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>           On entry, LDB specifies the first dimension of B as declared
+*>           in  the  calling  (sub)  program.   LDB  must  be  at  least
+*>           max( 1, m ).
+*> \endverbatim
+*>
+*> \param[in] BETA
+*> \verbatim
+*>          BETA is DOUBLE PRECISION.
+*>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
+*>           supplied as zero then C need not be set on input.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is DOUBLE PRECISION array, dimension ( LDC, N )
+*>           Before entry, the leading  m by n  part of the array  C must
+*>           contain the matrix  C,  except when  beta  is zero, in which
+*>           case C need not be set on entry.
+*>           On exit, the array  C  is overwritten by the  m by n updated
+*>           matrix.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>           On entry, LDC specifies the first dimension of C as declared
+*>           in  the  calling  (sub)  program.   LDC  must  be  at  least
+*>           max( 1, m ).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup skewhemm
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Level 3 Blas routine.
+*>  Derived from subroutine dsymm.
+*>
+*>  -- Written on 6-Jul-2025.
+*>     Shuo Zheng, China.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DSKEWSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,
+     +                     LDB,BETA,C,LDC)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level3 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION ALPHA,BETA
+      INTEGER LDA,LDB,LDC,M,N
+      CHARACTER SIDE,UPLO
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. External Functions ..
+      LOGICAL LSAME
+      EXTERNAL LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION TEMP1,TEMP2
+      INTEGER I,INFO,J,K,NROWA
+      LOGICAL UPPER
+*     ..
+*     .. Parameters ..
+      DOUBLE PRECISION ONE,ZERO
+      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+*     ..
+*
+*     Set NROWA as the number of rows of A.
+*
+      IF (LSAME(SIDE,'L')) THEN
+          NROWA = M
+      ELSE
+          NROWA = N
+      END IF
+      UPPER = LSAME(UPLO,'U')
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF ((.NOT.LSAME(SIDE,'L')) .AND.
+     +    (.NOT.LSAME(SIDE,'R'))) THEN
+          INFO = 1
+      ELSE IF ((.NOT.UPPER) .AND.
+     +         (.NOT.LSAME(UPLO,'L'))) THEN
+          INFO = 2
+      ELSE IF (M.LT.0) THEN
+          INFO = 3
+      ELSE IF (N.LT.0) THEN
+          INFO = 4
+      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
+          INFO = 7
+      ELSE IF (LDB.LT.MAX(1,M)) THEN
+          INFO = 9
+      ELSE IF (LDC.LT.MAX(1,M)) THEN
+          INFO = 12
+      END IF
+      IF (INFO.NE.0) THEN
+          CALL XERBLA('DSKEWSYMM ',INFO)
+          RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
+     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
+*
+*     And when  alpha.eq.zero.
+*
+      IF (ALPHA.EQ.ZERO) THEN
+          IF (BETA.EQ.ZERO) THEN
+              DO 20 J = 1,N
+                  DO 10 I = 1,M
+                      C(I,J) = ZERO
+   10             CONTINUE
+   20         CONTINUE
+          ELSE
+              DO 40 J = 1,N
+                  DO 30 I = 1,M
+                      C(I,J) = BETA*C(I,J)
+   30             CONTINUE
+   40         CONTINUE
+          END IF
+          RETURN
+      END IF
+*
+*     Start the operations.
+*
+      IF (LSAME(SIDE,'L')) THEN
+*
+*        Form  C := alpha*A*B + beta*C.
+*
+          IF (UPPER) THEN
+              DO 70 J = 1,N
+                  DO 60 I = 1,M
+                      TEMP1 = ALPHA*B(I,J)
+                      TEMP2 = ZERO
+                      DO 50 K = 1,I - 1
+                          C(K,J) = C(K,J) + TEMP1*A(K,I)
+                          TEMP2 = TEMP2 - B(K,J)*A(K,I)
+   50                 CONTINUE
+                      IF (BETA.EQ.ZERO) THEN
+                          C(I,J) = ALPHA*TEMP2
+                      ELSE
+                          C(I,J) = BETA*C(I,J) +
+     +                             ALPHA*TEMP2
+                      END IF
+   60             CONTINUE
+   70         CONTINUE
+          ELSE
+              DO 100 J = 1,N
+                  DO 90 I = M,1,-1
+                      TEMP1 = ALPHA*B(I,J)
+                      TEMP2 = ZERO
+                      DO 80 K = I + 1,M
+                          C(K,J) = C(K,J) + TEMP1*A(K,I)
+                          TEMP2 = TEMP2 - B(K,J)*A(K,I)
+   80                 CONTINUE
+                      IF (BETA.EQ.ZERO) THEN
+                          C(I,J) = ALPHA*TEMP2
+                      ELSE
+                          C(I,J) = BETA*C(I,J) +
+     +                             ALPHA*TEMP2
+                      END IF
+   90             CONTINUE
+  100         CONTINUE
+          END IF
+      ELSE
+*
+*        Form  C := alpha*B*A + beta*C.
+*
+          DO 170 J = 1,N
+              IF (BETA.EQ.ZERO) THEN
+                  DO 110 I = 1,M
+                      C(I,J) = ZERO
+  110             CONTINUE
+              ELSE
+                  DO 120 I = 1,M
+                      C(I,J) = BETA*C(I,J)
+  120             CONTINUE
+              END IF
+              DO 140 K = 1,J - 1
+                  IF (UPPER) THEN
+                      TEMP1 = ALPHA*A(K,J)
+                  ELSE
+                      TEMP1 = -ALPHA*A(J,K)
+                  END IF
+                  DO 130 I = 1,M
+                      C(I,J) = C(I,J) + TEMP1*B(I,K)
+  130             CONTINUE
+  140         CONTINUE
+              DO 160 K = J + 1,N
+                  IF (UPPER) THEN
+                      TEMP1 = -ALPHA*A(J,K)
+                  ELSE
+                      TEMP1 = ALPHA*A(K,J)
+                  END IF
+                  DO 150 I = 1,M
+                      C(I,J) = C(I,J) + TEMP1*B(I,K)
+  150             CONTINUE
+  160         CONTINUE
+  170     CONTINUE
+      END IF
+*
+      RETURN
+*
+*     End of DSKEWSYMM
+*
+      END
@@ -0,0 +1,327 @@
+*> \brief \b DSKEWSYMV
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DSKEWSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION ALPHA,BETA
+*       INTEGER INCX,INCY,LDA,N
+*       CHARACTER UPLO
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DSKEWSYMV performs the matrix-vector operation
+*>
+*>    y := alpha*A*x + beta*y,
+*>
+*> where alpha and beta are scalars, x and y are n element vectors and
+*> A is an n by n skew-symmetric matrix.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>           On entry, UPLO specifies whether the upper or lower
+*>           triangular part of the array A is to be referenced as
+*>           follows:
+*>
+*>              UPLO = 'U' or 'u'   Only the upper triangular part of A
+*>                                  is to be referenced.
+*>
+*>              UPLO = 'L' or 'l'   Only the lower triangular part of A
+*>                                  is to be referenced.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>           On entry, N specifies the order of the matrix A.
+*>           N must be at least zero.
+*> \endverbatim
+*>
+*> \param[in] ALPHA
+*> \verbatim
+*>          ALPHA is DOUBLE PRECISION
+*>           On entry, ALPHA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension ( LDA, N )
+*>           Before entry with  UPLO = 'U' or 'u', the strictly n by n
+*>           upper triangular part of the array A must contain the upper
+*>           triangular part of the skew-symmetric matrix and the leading
+*>           lower triangular part of A is not referenced.
+*>           Before entry with UPLO = 'L' or 'l', the strictly n by n
+*>           lower triangular part of the array A must contain the lower
+*>           triangular part of the skew-symmetric matrix and the leading
+*>           upper triangular part of A is not referenced.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>           On entry, LDA specifies the first dimension of A as declared
+*>           in the calling (sub) program. LDA must be at least
+*>           max( 1, n ).
+*> \endverbatim
+*>
+*> \param[in] X
+*> \verbatim
+*>          X is DOUBLE PRECISION array, dimension at least
+*>           ( 1 + ( n - 1 )*abs( INCX ) ).
+*>           Before entry, the incremented array X must contain the n
+*>           element vector x.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>           On entry, INCX specifies the increment for the elements of
+*>           X. INCX must not be zero.
+*> \endverbatim
+*>
+*> \param[in] BETA
+*> \verbatim
+*>          BETA is DOUBLE PRECISION.
+*>           On entry, BETA specifies the scalar beta. When BETA is
+*>           supplied as zero then Y need not be set on input.
+*> \endverbatim
+*>
+*> \param[in,out] Y
+*> \verbatim
+*>          Y is DOUBLE PRECISION array, dimension at least
+*>           ( 1 + ( n - 1 )*abs( INCY ) ).
+*>           Before entry, the incremented array Y must contain the n
+*>           element vector y. On exit, Y is overwritten by the updated
+*>           vector y.
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>           On entry, INCY specifies the increment for the elements of
+*>           Y. INCY must not be zero.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup skewhemv
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Level 2 Blas routine.
+*>  The vector and matrix arguments are not referenced when N = 0, or M = 0
+*>  Derived from subroutine dsymv.
+*>
+*>  -- Written on 6-Jul-2025.
+*>     Shuo Zheng, China.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DSKEWSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level2 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION ALPHA,BETA
+      INTEGER INCX,INCY,LDA,N
+      CHARACTER UPLO
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ONE,ZERO
+      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION TEMP1,TEMP2
+      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
+*     ..
+*     .. External Functions ..
+      LOGICAL LSAME
+      EXTERNAL LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX
+*     ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
+          INFO = 1
+      ELSE IF (N.LT.0) THEN
+          INFO = 2
+      ELSE IF (LDA.LT.MAX(1,N)) THEN
+          INFO = 5
+      ELSE IF (INCX.EQ.0) THEN
+          INFO = 7
+      ELSE IF (INCY.EQ.0) THEN
+          INFO = 10
+      END IF
+      IF (INFO.NE.0) THEN
+          CALL XERBLA('DSKEWSYMV ',INFO)
+          RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
+*
+*     Set up the start points in  X  and  Y.
+*
+      IF (INCX.GT.0) THEN
+          KX = 1
+      ELSE
+          KX = 1 - (N-1)*INCX
+      END IF
+      IF (INCY.GT.0) THEN
+          KY = 1
+      ELSE
+          KY = 1 - (N-1)*INCY
+      END IF
+*
+*     Start the operations. In this version the elements of A are
+*     accessed sequentially with one pass through the triangular part
+*     of A.
+*
+*     First form  y := beta*y.
+*
+      IF (BETA.NE.ONE) THEN
+          IF (INCY.EQ.1) THEN
+              IF (BETA.EQ.ZERO) THEN
+                  DO 10 I = 1,N
+                      Y(I) = ZERO
+   10             CONTINUE
+              ELSE
+                  DO 20 I = 1,N
+                      Y(I) = BETA*Y(I)
+   20             CONTINUE
+              END IF
+          ELSE
+              IY = KY
+              IF (BETA.EQ.ZERO) THEN
+                  DO 30 I = 1,N
+                      Y(IY) = ZERO
+                      IY = IY + INCY
+   30             CONTINUE
+              ELSE
+                  DO 40 I = 1,N
+                      Y(IY) = BETA*Y(IY)
+                      IY = IY + INCY
+   40             CONTINUE
+              END IF
+          END IF
+      END IF
+      IF (ALPHA.EQ.ZERO) RETURN
+      IF (LSAME(UPLO,'U')) THEN
+*
+*        Form  y  when A is stored in upper triangle.
+*
+          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
+              DO 60 J = 1,N
+                  TEMP1 = ALPHA*X(J)
+                  TEMP2 = ZERO
+                  DO 50 I = 1,J - 1
+                      Y(I) = Y(I) + TEMP1*A(I,J)
+                      TEMP2 = TEMP2 - A(I,J)*X(I)
+   50             CONTINUE
+                  Y(J) = Y(J) + ALPHA*TEMP2
+   60         CONTINUE
+          ELSE
+              JX = KX
+              JY = KY
+              DO 80 J = 1,N
+                  TEMP1 = ALPHA*X(JX)
+                  TEMP2 = ZERO
+                  IX = KX
+                  IY = KY
+                  DO 70 I = 1,J - 1
+                      Y(IY) = Y(IY) + TEMP1*A(I,J)
+                      TEMP2 = TEMP2 - A(I,J)*X(IX)
+                      IX = IX + INCX
+                      IY = IY + INCY
+   70             CONTINUE
+                  Y(JY) = Y(JY) + ALPHA*TEMP2
+                  JX = JX + INCX
+                  JY = JY + INCY
+   80         CONTINUE
+          END IF
+      ELSE
+*
+*        Form  y  when A is stored in lower triangle.
+*
+          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
+              DO 100 J = 1,N
+                  TEMP1 = ALPHA*X(J)
+                  TEMP2 = ZERO
+                  DO 90 I = J + 1,N
+                      Y(I) = Y(I) + TEMP1*A(I,J)
+                      TEMP2 = TEMP2 - A(I,J)*X(I)
+   90             CONTINUE
+                  Y(J) = Y(J) + ALPHA*TEMP2
+  100         CONTINUE
+          ELSE
+              JX = KX
+              JY = KY
+              DO 120 J = 1,N
+                  TEMP1 = ALPHA*X(JX)
+                  TEMP2 = ZERO
+                  IX = JX
+                  IY = JY
+                  DO 110 I = J + 1,N
+                      IX = IX + INCX
+                      IY = IY + INCY
+                      Y(IY) = Y(IY) + TEMP1*A(I,J)
+                      TEMP2 = TEMP2 - A(I,J)*X(IX)
+  110             CONTINUE
+                  Y(JY) = Y(JY) + ALPHA*TEMP2
+                  JX = JX + INCX
+                  JY = JY + INCY
+  120         CONTINUE
+          END IF
+      END IF
+*
+      RETURN
+*
+*     End of DSKEWSYMV
+*
+      END
@@ -0,0 +1,294 @@
+*> \brief \b DSKEWSYR2
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DSKEWSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION ALPHA
+*       INTEGER INCX,INCY,LDA,N
+*       CHARACTER UPLO
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DSKEWSYR2 performs the skew-symmetric rank 2 operation
+*>
+*>    A := -alpha*x*y**T + alpha*y*x**T + A,
+*>
+*> where alpha is a scalar, x and y are n element vectors and A is an n
+*> by n skew-symmetric matrix.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>           On entry, UPLO specifies whether the upper or lower
+*>           triangular part of the array A is to be referenced as
+*>           follows:
+*>
+*>              UPLO = 'U' or 'u'   Only the upper triangular part of A
+*>                                  is to be referenced.
+*>
+*>              UPLO = 'L' or 'l'   Only the lower triangular part of A
+*>                                  is to be referenced.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>           On entry, N specifies the order of the matrix A.
+*>           N must be at least zero.
+*> \endverbatim
+*>
+*> \param[in] ALPHA
+*> \verbatim
+*>          ALPHA is DOUBLE PRECISION
+*>           On entry, ALPHA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] X
+*> \verbatim
+*>          X is DOUBLE PRECISION array, dimension at least
+*>           ( 1 + ( n - 1 )*abs( INCX ) ).
+*>           Before entry, the incremented array X must contain the n
+*>           element vector x.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>           On entry, INCX specifies the increment for the elements of
+*>           X. INCX must not be zero.
+*> \endverbatim
+*>
+*> \param[in] Y
+*> \verbatim
+*>          Y is DOUBLE PRECISION array, dimension at least
+*>           ( 1 + ( n - 1 )*abs( INCY ) ).
+*>           Before entry, the incremented array Y must contain the n
+*>           element vector y.
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>           On entry, INCY specifies the increment for the elements of
+*>           Y. INCY must not be zero.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension ( LDA, N )
+*>           Before entry with  UPLO = 'U' or 'u', the strictly n by n
+*>           upper triangular part of the array A must contain the upper
+*>           triangular part of the skew-symmetric matrix and the leading
+*>           lower triangular part of A is not referenced. On exit, the
+*>           upper triangular part of the array A is overwritten by the
+*>           upper triangular part of the updated matrix.
+*>           Before entry with UPLO = 'L' or 'l', the strictly n by n
+*>           lower triangular part of the array A must contain the lower
+*>           triangular part of the skew-symmetric matrix and the leading
+*>           upper triangular part of A is not referenced. On exit, the
+*>           lower triangular part of the array A is overwritten by the
+*>           lower triangular part of the updated matrix.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>           On entry, LDA specifies the first dimension of A as declared
+*>           in the calling (sub) program. LDA must be at least
+*>           max( 1, n ).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup skewher2
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Level 2 Blas routine.
+*>  Derived from subroutine dsyr2.
+*>
+*>  -- Written on 6-Jul-2025.
+*>     Shuo Zheng, China.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DSKEWSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level2 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION ALPHA
+      INTEGER INCX,INCY,LDA,N
+      CHARACTER UPLO
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION TEMP1,TEMP2
+      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
+*     ..
+*     .. External Functions ..
+      LOGICAL LSAME
+      EXTERNAL LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX
+*     ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
+          INFO = 1
+      ELSE IF (N.LT.0) THEN
+          INFO = 2
+      ELSE IF (INCX.EQ.0) THEN
+          INFO = 5
+      ELSE IF (INCY.EQ.0) THEN
+          INFO = 7
+      ELSE IF (LDA.LT.MAX(1,N)) THEN
+          INFO = 9
+      END IF
+      IF (INFO.NE.0) THEN
+          CALL XERBLA('DSKEWSYR2 ',INFO)
+          RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
+*
+*     Set up the start points in X and Y if the increments are not both
+*     unity.
+*
+      IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
+          IF (INCX.GT.0) THEN
+              KX = 1
+          ELSE
+              KX = 1 - (N-1)*INCX
+          END IF
+          IF (INCY.GT.0) THEN
+              KY = 1
+          ELSE
+              KY = 1 - (N-1)*INCY
+          END IF
+          JX = KX
+          JY = KY
+      END IF
+*
+*     Start the operations. In this version the elements of A are
+*     accessed sequentially with one pass through the triangular part
+*     of A.
+*
+      IF (LSAME(UPLO,'U')) THEN
+*
+*        Form  A  when A is stored in the upper triangle.
+*
+          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
+              DO 20 J = 1,N
+                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
+                      TEMP1 = ALPHA*Y(J)
+                      TEMP2 = ALPHA*X(J)
+                      DO 10 I = 1,J-1
+                          A(I,J) = A(I,J) - X(I)*TEMP1 + Y(I)*TEMP2
+   10                 CONTINUE
+                  END IF
+   20         CONTINUE
+          ELSE
+              DO 40 J = 1,N
+                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
+                      TEMP1 = ALPHA*Y(JY)
+                      TEMP2 = ALPHA*X(JX)
+                      IX = KX
+                      IY = KY
+                      DO 30 I = 1,J-1
+                          A(I,J) = A(I,J) - X(IX)*TEMP1 + Y(IY)*TEMP2
+                          IX = IX + INCX
+                          IY = IY + INCY
+   30                 CONTINUE
+                  END IF
+                  JX = JX + INCX
+                  JY = JY + INCY
+   40         CONTINUE
+          END IF
+      ELSE
+*
+*        Form  A  when A is stored in the lower triangle.
+*
+          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
+              DO 60 J = 1,N
+                  IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
+                      TEMP1 = ALPHA*Y(J)
+                      TEMP2 = ALPHA*X(J)
+                      DO 50 I = J+1,N
+                          A(I,J) = A(I,J) - X(I)*TEMP1 + Y(I)*TEMP2
+   50                 CONTINUE
+                  END IF
+   60         CONTINUE
+          ELSE
+              DO 80 J = 1,N
+                  IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
+                      TEMP1 = ALPHA*Y(JY)
+                      TEMP2 = ALPHA*X(JX)
+                      IX = JX + INCX
+                      IY = JY + INCY
+                      DO 70 I = J+1,N
+                          A(I,J) = A(I,J) - X(IX)*TEMP1 + Y(IY)*TEMP2
+                          IX = IX + INCX
+                          IY = IY + INCY
+   70                 CONTINUE
+                  END IF
+                  JX = JX + INCX
+                  JY = JY + INCY
+   80         CONTINUE
+          END IF
+      END IF
+*
+      RETURN
+*
+*     End of DSKEWSYR2
+*
+      END
@@ -0,0 +1,395 @@
+*> \brief \b DSKEWSYR2K
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DSKEWSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+*
+*       .. Scalar Arguments ..
+*       DOUBLE PRECISION ALPHA,BETA
+*       INTEGER K,LDA,LDB,LDC,N
+*       CHARACTER TRANS,UPLO
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DSKEWSYR2K performs one of the skew-symmetric rank 2k operations
+*>
+*>    C := -alpha*A*B**T + alpha*B*A**T + beta*C,
+*>
+*> or
+*>
+*>    C := -alpha*A**T*B + alpha*B**T*A + beta*C,
+*>
+*> where  alpha and beta  are scalars, C is an  n by n  skew-symmetric matrix
+*> and  A and B  are  n by k  matrices  in the  first  case  and  k by n
+*> matrices in the second case.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
+*>           triangular  part  of the  array  C  is to be  referenced  as
+*>           follows:
+*>
+*>              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
+*>                                  is to be referenced.
+*>
+*>              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
+*>                                  is to be referenced.
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>           On entry,  TRANS  specifies the operation to be performed as
+*>           follows:
+*>
+*>              TRANS = 'N' or 'n'   C := -alpha*A*B**T + alpha*B*A**T +
+*>                                        beta*C.
+*>
+*>              TRANS = 'T' or 't'   C := -alpha*A**T*B + alpha*B**T*A +
+*>                                        beta*C.
+*>
+*>              TRANS = 'C' or 'c'   C := -alpha*A**T*B + alpha*B**T*A +
+*>                                        beta*C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>           On entry,  N specifies the order of the matrix C.  N must be
+*>           at least zero.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
+*>           of  columns  of the  matrices  A and B,  and on  entry  with
+*>           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
+*>           of rows of the matrices  A and B.  K must be at least  zero.
+*> \endverbatim
+*>
+*> \param[in] ALPHA
+*> \verbatim
+*>          ALPHA is DOUBLE PRECISION.
+*>           On entry, ALPHA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*>          A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
+*>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
+*>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
+*>           part of the array  A  must contain the matrix  A,  otherwise
+*>           the leading  k by n  part of the array  A  must contain  the
+*>           matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>           On entry, LDA specifies the first dimension of A as declared
+*>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
+*>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
+*>           be at least  max( 1, k ).
+*> \endverbatim
+*>
+*> \param[in] B
+*> \verbatim
+*>          B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
+*>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
+*>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
+*>           part of the array  B  must contain the matrix  B,  otherwise
+*>           the leading  k by n  part of the array  B  must contain  the
+*>           matrix B.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*>          LDB is INTEGER
+*>           On entry, LDB specifies the first dimension of B as declared
+*>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
+*>           then  LDB must be at least  max( 1, n ), otherwise  LDB must
+*>           be at least  max( 1, k ).
+*> \endverbatim
+*>
+*> \param[in] BETA
+*> \verbatim
+*>          BETA is DOUBLE PRECISION.
+*>           On entry, BETA specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is DOUBLE PRECISION array, dimension ( LDC, N )
+*>           Before entry  with  UPLO = 'U' or 'u',  the strictly  n by n
+*>           upper triangular part of the array C must contain the upper
+*>           triangular part  of the  skew-symmetric matrix  and the leading
+*>           lower triangular part of C is not referenced.  On exit, the
+*>           upper triangular part of the array  C is overwritten by the
+*>           upper triangular part of the updated matrix.
+*>           Before entry  with  UPLO = 'L' or 'l',  the strictly  n by n
+*>           lower triangular part of the array C must contain the lower
+*>           triangular part  of the  skew-symmetric matrix  and the leading
+*>           upper triangular part of C is not referenced.  On exit, the
+*>           lower triangular part of the array  C is overwritten by the
+*>           lower triangular part of the updated matrix.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>           On entry, LDC specifies the first dimension of C as declared
+*>           in  the  calling  (sub)  program.   LDC  must  be  at  least
+*>           max( 1, n ).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup skewher2k
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  Level 3 Blas routine.
+*>  Derived from subroutine dsyr2k.
+*>
+*>  -- Written on 6-Jul-2025.
+*>     Shuo Zheng, China.
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DSKEWSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,
+     +                      LDB,BETA,C,LDC)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level3 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION ALPHA,BETA
+      INTEGER K,LDA,LDB,LDC,N
+      CHARACTER TRANS,UPLO
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+*     ..
+*
+*  =====================================================================
+*
+*     .. External Functions ..
+      LOGICAL LSAME
+      EXTERNAL LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MAX
+*     ..
+*     .. Local Scalars ..
+      DOUBLE PRECISION TEMP1,TEMP2
+      INTEGER I,INFO,J,L,NROWA
+      LOGICAL UPPER
+*     ..
+*     .. Parameters ..
+      DOUBLE PRECISION ONE,ZERO
+      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+*     ..
+*
+*     Test the input parameters.
+*
+      IF (LSAME(TRANS,'N')) THEN
+          NROWA = N
+      ELSE
+          NROWA = K
+      END IF
+      UPPER = LSAME(UPLO,'U')
+*
+      INFO = 0
+      IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
+          INFO = 1
+      ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
+     +         (.NOT.LSAME(TRANS,'T')) .AND.
+     +         (.NOT.LSAME(TRANS,'C'))) THEN
+          INFO = 2
+      ELSE IF (N.LT.0) THEN
+          INFO = 3
+      ELSE IF (K.LT.0) THEN
+          INFO = 4
+      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
+          INFO = 7
+      ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
+          INFO = 9
+      ELSE IF (LDC.LT.MAX(1,N)) THEN
+          INFO = 12
+      END IF
+      IF (INFO.NE.0) THEN
+          CALL XERBLA('DSKEWSYR2K',INFO)
+          RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
+     +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
+*
+*     And when  alpha.eq.zero.
+*
+      IF (ALPHA.EQ.ZERO) THEN
+          IF (UPPER) THEN
+              IF (BETA.EQ.ZERO) THEN
+                  DO 20 J = 1,N
+                      DO 10 I = 1,J-1
+                          C(I,J) = ZERO
+   10                 CONTINUE
+   20             CONTINUE
+              ELSE
+                  DO 40 J = 1,N
+                      DO 30 I = 1,J-1
+                          C(I,J) = BETA*C(I,J)
+   30                 CONTINUE
+   40             CONTINUE
+              END IF
+          ELSE
+              IF (BETA.EQ.ZERO) THEN
+                  DO 60 J = 1,N
+                      DO 50 I = J+1,N
+                          C(I,J) = ZERO
+   50                 CONTINUE
+   60             CONTINUE
+              ELSE
+                  DO 80 J = 1,N
+                      DO 70 I = J+1,N
+                          C(I,J) = BETA*C(I,J)
+   70                 CONTINUE
+   80             CONTINUE
+              END IF
+          END IF
+          RETURN
+      END IF
+*
+*     Start the operations.
+*
+      IF (LSAME(TRANS,'N')) THEN
+*
+*        Form  C := alpha*A*B**T + alpha*B*A**T + C.
+*
+          IF (UPPER) THEN
+              DO 130 J = 1,N
+                  IF (BETA.EQ.ZERO) THEN
+                      DO 90 I = 1,J-1
+                          C(I,J) = ZERO
+   90                 CONTINUE
+                  ELSE IF (BETA.NE.ONE) THEN
+                      DO 100 I = 1,J-1
+                          C(I,J) = BETA*C(I,J)
+  100                 CONTINUE
+                  END IF
+                  DO 120 L = 1,K
+                      IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
+                          TEMP1 = ALPHA*B(J,L)
+                          TEMP2 = ALPHA*A(J,L)
+                          DO 110 I = 1,J-1
+                              C(I,J) = C(I,J) - A(I,L)*TEMP1 +
+     +                                 B(I,L)*TEMP2
+  110                     CONTINUE
+                      END IF
+  120             CONTINUE
+  130         CONTINUE
+          ELSE
+              DO 180 J = 1,N
+                  IF (BETA.EQ.ZERO) THEN
+                      DO 140 I = J+1,N
+                          C(I,J) = ZERO
+  140                 CONTINUE
+                  ELSE IF (BETA.NE.ONE) THEN
+                      DO 150 I = J+1,N
+                          C(I,J) = BETA*C(I,J)
+  150                 CONTINUE
+                  END IF
+                  DO 170 L = 1,K
+                      IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
+                          TEMP1 = ALPHA*B(J,L)
+                          TEMP2 = ALPHA*A(J,L)
+                          DO 160 I = J+1,N
+                              C(I,J) = C(I,J) - A(I,L)*TEMP1 +
+     +                                 B(I,L)*TEMP2
+  160                     CONTINUE
+                      END IF
+  170             CONTINUE
+  180         CONTINUE
+          END IF
+      ELSE
+*
+*        Form  C := alpha*A**T*B + alpha*B**T*A + C.
+*
+          IF (UPPER) THEN
+              DO 210 J = 1,N
+                  DO 200 I = 1,J-1
+                      TEMP1 = ZERO
+                      TEMP2 = ZERO
+                      DO 190 L = 1,K
+                          TEMP1 = TEMP1 + A(L,I)*B(L,J)
+                          TEMP2 = TEMP2 + B(L,I)*A(L,J)
+  190                 CONTINUE
+                      IF (BETA.EQ.ZERO) THEN
+                          C(I,J) = -ALPHA*TEMP1 + ALPHA*TEMP2
+                      ELSE
+                          C(I,J) = BETA*C(I,J) - ALPHA*TEMP1 +
+     +                             ALPHA*TEMP2
+                      END IF
+  200             CONTINUE
+  210         CONTINUE
+          ELSE
+              DO 240 J = 1,N
+                  DO 230 I = J+1,N
+                      TEMP1 = ZERO
+                      TEMP2 = ZERO
+                      DO 220 L = 1,K
+                          TEMP1 = TEMP1 + A(L,I)*B(L,J)
+                          TEMP2 = TEMP2 + B(L,I)*A(L,J)
+  220                 CONTINUE
+                      IF (BETA.EQ.ZERO) THEN
+                          C(I,J) = -ALPHA*TEMP1 + ALPHA*TEMP2
+                      ELSE
+                          C(I,J) = BETA*C(I,J) - ALPHA*TEMP1 +
+     +                             ALPHA*TEMP2
+                      END IF
+  230             CONTINUE
+  240         CONTINUE
+          END IF
+      END IF
+*
+      RETURN
+*
+*     End of DSKEWSYR2K
+*
+      END
@@ -144,6 +144,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -124,6 +124,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -139,6 +139,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -79,6 +79,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSWAP(N,DX,INCX,DY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -186,6 +186,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -149,6 +149,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -129,6 +129,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYR(UPLO,N,ALPHA,X,INCX,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -144,6 +144,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -189,6 +189,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -166,6 +166,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -183,6 +183,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -197,10 +198,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -270,28 +267,24 @@
              KPLUS1 = K + 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          L = KPLUS1 - J
-                          DO 10 I = MAX(1,J-K),J - 1
-                              X(I) = X(I) + TEMP*A(L+I,J)
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
-                      END IF
+                      TEMP = X(J)
+                      L = KPLUS1 - J
+                      DO 10 I = MAX(1,J-K),J - 1
+                          X(I) = X(I) + TEMP*A(L+I,J)
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          L = KPLUS1 - J
-                          DO 30 I = MAX(1,J-K),J - 1
-                              X(IX) = X(IX) + TEMP*A(L+I,J)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      L = KPLUS1 - J
+                      DO 30 I = MAX(1,J-K),J - 1
+                          X(IX) = X(IX) + TEMP*A(L+I,J)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
                      JX = JX + INCX
                      IF (J.GT.K) KX = KX + INCX
   40             CONTINUE
@@ -299,29 +292,25 @@
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          L = 1 - J
-                          DO 50 I = MIN(N,J+K),J + 1,-1
-                              X(I) = X(I) + TEMP*A(L+I,J)
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(1,J)
-                      END IF
+                      TEMP = X(J)
+                      L = 1 - J
+                      DO 50 I = MIN(N,J+K),J + 1,-1
+                          X(I) = X(I) + TEMP*A(L+I,J)
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(1,J)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          L = 1 - J
-                          DO 70 I = MIN(N,J+K),J + 1,-1
-                              X(IX) = X(IX) + TEMP*A(L+I,J)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(1,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      L = 1 - J
+                      DO 70 I = MIN(N,J+K),J + 1,-1
+                          X(IX) = X(IX) + TEMP*A(L+I,J)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(1,J)
                      JX = JX - INCX
                      IF ((N-J).GE.K) KX = KX - INCX
   80             CONTINUE
@@ -186,6 +186,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -200,10 +201,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -273,59 +270,51 @@
              KPLUS1 = K + 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          L = KPLUS1 - J
-                          IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
-                          TEMP = X(J)
-                          DO 10 I = J - 1,MAX(1,J-K),-1
-                              X(I) = X(I) - TEMP*A(L+I,J)
-   10                     CONTINUE
-                      END IF
+                      L = KPLUS1 - J
+                      IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
+                      TEMP = X(J)
+                      DO 10 I = J - 1,MAX(1,J-K),-1
+                          X(I) = X(I) - TEMP*A(L+I,J)
+   10                 CONTINUE
   20             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 40 J = N,1,-1
                      KX = KX - INCX
-                      IF (X(JX).NE.ZERO) THEN
-                          IX = KX
-                          L = KPLUS1 - J
-                          IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
-                          TEMP = X(JX)
-                          DO 30 I = J - 1,MAX(1,J-K),-1
-                              X(IX) = X(IX) - TEMP*A(L+I,J)
-                              IX = IX - INCX
-   30                     CONTINUE
-                      END IF
+                      IX = KX
+                      L = KPLUS1 - J
+                      IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
+                      TEMP = X(JX)
+                      DO 30 I = J - 1,MAX(1,J-K),-1
+                          X(IX) = X(IX) - TEMP*A(L+I,J)
+                          IX = IX - INCX
+   30                 CONTINUE
                      JX = JX - INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          L = 1 - J
-                          IF (NOUNIT) X(J) = X(J)/A(1,J)
-                          TEMP = X(J)
-                          DO 50 I = J + 1,MIN(N,J+K)
-                              X(I) = X(I) - TEMP*A(L+I,J)
-   50                     CONTINUE
-                      END IF
+                      L = 1 - J
+                      IF (NOUNIT) X(J) = X(J)/A(1,J)
+                      TEMP = X(J)
+                      DO 50 I = J + 1,MIN(N,J+K)
+                          X(I) = X(I) - TEMP*A(L+I,J)
+   50                 CONTINUE
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
                      KX = KX + INCX
-                      IF (X(JX).NE.ZERO) THEN
-                          IX = KX
-                          L = 1 - J
-                          IF (NOUNIT) X(JX) = X(JX)/A(1,J)
-                          TEMP = X(JX)
-                          DO 70 I = J + 1,MIN(N,J+K)
-                              X(IX) = X(IX) - TEMP*A(L+I,J)
-                              IX = IX + INCX
-   70                     CONTINUE
-                      END IF
+                      IX = KX
+                      L = 1 - J
+                      IF (NOUNIT) X(JX) = X(JX)/A(1,J)
+                      TEMP = X(JX)
+                      DO 70 I = J + 1,MIN(N,J+K)
+                          X(IX) = X(IX) - TEMP*A(L+I,J)
+                          IX = IX + INCX
+   70                 CONTINUE
                      JX = JX + INCX
   80             CONTINUE
              END IF
@@ -139,6 +139,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -153,10 +154,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -219,29 +216,25 @@
              KK = 1
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          K = KK
-                          DO 10 I = 1,J - 1
-                              X(I) = X(I) + TEMP*AP(K)
-                              K = K + 1
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
-                      END IF
+                      TEMP = X(J)
+                      K = KK
+                      DO 10 I = 1,J - 1
+                          X(I) = X(I) + TEMP*AP(K)
+                          K = K + 1
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
                      KK = KK + J
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 30 K = KK,KK + J - 2
-                              X(IX) = X(IX) + TEMP*AP(K)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 30 K = KK,KK + J - 2
+                          X(IX) = X(IX) + TEMP*AP(K)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
                      JX = JX + INCX
                      KK = KK + J
   40             CONTINUE
@@ -250,30 +243,26 @@
              KK = (N* (N+1))/2
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          K = KK
-                          DO 50 I = N,J + 1,-1
-                              X(I) = X(I) + TEMP*AP(K)
-                              K = K - 1
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
-                      END IF
+                      TEMP = X(J)
+                      K = KK
+                      DO 50 I = N,J + 1,-1
+                          X(I) = X(I) + TEMP*AP(K)
+                          K = K - 1
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
                      KK = KK - (N-J+1)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 70 K = KK,KK - (N- (J+1)),-1
-                              X(IX) = X(IX) + TEMP*AP(K)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 70 K = KK,KK - (N- (J+1)),-1
+                          X(IX) = X(IX) + TEMP*AP(K)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
                      JX = JX - INCX
                      KK = KK - (N-J+1)
   80             CONTINUE
@@ -141,6 +141,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -155,10 +156,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -221,29 +218,25 @@
              KK = (N* (N+1))/2
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK - 1
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K - 1
-   10                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/AP(KK)
+                      TEMP = X(J)
+                      K = KK - 1
+                      DO 10 I = J - 1,1,-1
+                          X(I) = X(I) - TEMP*AP(K)
+                          K = K - 1
+   10                 CONTINUE
                      KK = KK - J
   20             CONTINUE
              ELSE
                  JX = KX + (N-1)*INCX
                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 K = KK - 1,KK - J + 1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   30                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 30 K = KK - 1,KK - J + 1,-1
+                          IX = IX - INCX
+                          X(IX) = X(IX) - TEMP*AP(K)
+   30                 CONTINUE
                      JX = JX - INCX
                      KK = KK - J
   40             CONTINUE
@@ -252,29 +245,25 @@
              KK = 1
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/AP(KK)
-                          TEMP = X(J)
-                          K = KK + 1
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*AP(K)
-                              K = K + 1
-   50                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/AP(KK)
+                      TEMP = X(J)
+                      K = KK + 1
+                      DO 50 I = J + 1,N
+                          X(I) = X(I) - TEMP*AP(K)
+                          K = K + 1
+   50                 CONTINUE
                      KK = KK + (N-J+1)
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/AP(KK)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 K = KK + 1,KK + N - J
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*AP(K)
-   70                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 70 K = KK + 1,KK + N - J
+                          IX = IX + INCX
+                          X(IX) = X(IX) - TEMP*AP(K)
+   70                 CONTINUE
                      JX = JX + INCX
                      KK = KK + (N-J+1)
   80             CONTINUE
@@ -174,6 +174,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -272,27 +273,23 @@
              IF (UPPER) THEN
                  DO 50 J = 1,N
                      DO 40 K = 1,M
-                          IF (B(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*B(K,J)
-                              DO 30 I = 1,K - 1
-                                  B(I,J) = B(I,J) + TEMP*A(I,K)
-   30                         CONTINUE
-                              IF (NOUNIT) TEMP = TEMP*A(K,K)
-                              B(K,J) = TEMP
-                          END IF
+                          TEMP = ALPHA*B(K,J)
+                          DO 30 I = 1,K - 1
+                              B(I,J) = B(I,J) + TEMP*A(I,K)
+   30                     CONTINUE
+                          IF (NOUNIT) TEMP = TEMP*A(K,K)
+                          B(K,J) = TEMP
   40                 CONTINUE
   50             CONTINUE
              ELSE
                  DO 80 J = 1,N
                      DO 70 K = M,1,-1
-                          IF (B(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*B(K,J)
-                              B(K,J) = TEMP
-                              IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
-                              DO 60 I = K + 1,M
-                                  B(I,J) = B(I,J) + TEMP*A(I,K)
-   60                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*B(K,J)
+                          B(K,J) = TEMP
+                          IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
+                          DO 60 I = K + 1,M
+                              B(I,J) = B(I,J) + TEMP*A(I,K)
+   60                     CONTINUE
   70                 CONTINUE
   80             CONTINUE
              END IF
@@ -337,12 +334,10 @@
                          B(I,J) = TEMP*B(I,J)
  150                 CONTINUE
                      DO 170 K = 1,J - 1
-                          IF (A(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*A(K,J)
-                              DO 160 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  160                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(K,J)
+                          DO 160 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  160                     CONTINUE
  170                 CONTINUE
  180             CONTINUE
              ELSE
@@ -353,12 +348,10 @@
                          B(I,J) = TEMP*B(I,J)
  190                 CONTINUE
                      DO 210 K = J + 1,N
-                          IF (A(K,J).NE.ZERO) THEN
-                              TEMP = ALPHA*A(K,J)
-                              DO 200 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  200                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(K,J)
+                          DO 200 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  200                     CONTINUE
  210                 CONTINUE
  220             CONTINUE
              END IF
@@ -369,12 +362,10 @@
              IF (UPPER) THEN
                  DO 260 K = 1,N
                      DO 240 J = 1,K - 1
-                          IF (A(J,K).NE.ZERO) THEN
-                              TEMP = ALPHA*A(J,K)
-                              DO 230 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  230                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(J,K)
+                          DO 230 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  230                     CONTINUE
  240                 CONTINUE
                      TEMP = ALPHA
                      IF (NOUNIT) TEMP = TEMP*A(K,K)
@@ -387,12 +378,10 @@
              ELSE
                  DO 300 K = N,1,-1
                      DO 280 J = K + 1,N
-                          IF (A(J,K).NE.ZERO) THEN
-                              TEMP = ALPHA*A(J,K)
-                              DO 270 I = 1,M
-                                  B(I,J) = B(I,J) + TEMP*B(I,K)
-  270                         CONTINUE
-                          END IF
+                          TEMP = ALPHA*A(J,K)
+                          DO 270 I = 1,M
+                              B(I,J) = B(I,J) + TEMP*B(I,K)
+  270                     CONTINUE
  280                 CONTINUE
                      TEMP = ALPHA
                      IF (NOUNIT) TEMP = TEMP*A(K,K)
@@ -144,6 +144,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -158,10 +159,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -228,53 +225,45 @@
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          DO 10 I = 1,J - 1
-                              X(I) = X(I) + TEMP*A(I,J)
-   10                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(J,J)
-                      END IF
+                      TEMP = X(J)
+                      DO 10 I = 1,J - 1
+                          X(I) = X(I) + TEMP*A(I,J)
+   10                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(J,J)
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 30 I = 1,J - 1
-                              X(IX) = X(IX) + TEMP*A(I,J)
-                              IX = IX + INCX
-   30                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 30 I = 1,J - 1
+                          X(IX) = X(IX) + TEMP*A(I,J)
+                          IX = IX + INCX
+   30                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      JX = JX + INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          TEMP = X(J)
-                          DO 50 I = N,J + 1,-1
-                              X(I) = X(I) + TEMP*A(I,J)
-   50                     CONTINUE
-                          IF (NOUNIT) X(J) = X(J)*A(J,J)
-                      END IF
+                      TEMP = X(J)
+                      DO 50 I = N,J + 1,-1
+                          X(I) = X(I) + TEMP*A(I,J)
+   50                 CONTINUE
+                      IF (NOUNIT) X(J) = X(J)*A(J,J)
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          TEMP = X(JX)
-                          IX = KX
-                          DO 70 I = N,J + 1,-1
-                              X(IX) = X(IX) + TEMP*A(I,J)
-                              IX = IX - INCX
-   70                     CONTINUE
-                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
-                      END IF
+                      TEMP = X(JX)
+                      IX = KX
+                      DO 70 I = N,J + 1,-1
+                          X(IX) = X(IX) + TEMP*A(I,J)
+                          IX = IX - INCX
+   70                 CONTINUE
+                      IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      JX = JX - INCX
   80             CONTINUE
              END IF
@@ -178,6 +178,7 @@
 *>
 *  =====================================================================
      SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level3 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -210,8 +211,8 @@
      LOGICAL LSIDE,NOUNIT,UPPER
 *     ..
 *     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
-      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *
 *     Test the input parameters.
@@ -275,34 +276,26 @@
 *
              IF (UPPER) THEN
                  DO 60 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 30 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-   30                     CONTINUE
-                      END IF
-                      DO 50 K = M,1,-1
-                          IF (B(K,J).NE.ZERO) THEN
-                              IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
-                              DO 40 I = 1,K - 1
-                                  B(I,J) = B(I,J) - B(K,J)*A(I,K)
-   40                         CONTINUE
-                          END IF
+                       DO 30 I = 1,M
+                           B(I,J) = ALPHA*B(I,J)
+   30                  CONTINUE
+                       DO 50 K = M,1,-1
+                           IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+                           DO 40 I = 1,K - 1
+                               B(I,J) = B(I,J) - B(K,J)*A(I,K)
+   40                     CONTINUE
   50                 CONTINUE
   60             CONTINUE
              ELSE
                  DO 100 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 70 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-   70                     CONTINUE
-                      END IF
-                      DO 90 K = 1,M
-                          IF (B(K,J).NE.ZERO) THEN
-                              IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
-                              DO 80 I = K + 1,M
-                                  B(I,J) = B(I,J) - B(K,J)*A(I,K)
-   80                         CONTINUE
-                          END IF
+                       DO 70 I = 1,M
+                           B(I,J) = ALPHA*B(I,J)
+   70                  CONTINUE
+                       DO 90 K = 1,M
+                           IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+                           DO 80 I = K + 1,M
+                               B(I,J) = B(I,J) - B(K,J)*A(I,K)
+   80                     CONTINUE
   90                 CONTINUE
  100             CONTINUE
              END IF
@@ -341,43 +334,33 @@
 *
              IF (UPPER) THEN
                  DO 210 J = 1,N
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 170 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-  170                     CONTINUE
-                      END IF
+                      DO 170 I = 1,M
+                          B(I,J) = ALPHA*B(I,J)
+  170                 CONTINUE
                      DO 190 K = 1,J - 1
-                          IF (A(K,J).NE.ZERO) THEN
-                              DO 180 I = 1,M
-                                  B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  180                         CONTINUE
-                          END IF
+                          DO 180 I = 1,M
+                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
+  180                     CONTINUE
  190                 CONTINUE
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(J,J)
                          DO 200 I = 1,M
-                              B(I,J) = TEMP*B(I,J)
+                              B(I,J) = B(I,J)/A(J,J)
  200                     CONTINUE
                      END IF
  210             CONTINUE
              ELSE
                  DO 260 J = N,1,-1
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 220 I = 1,M
-                              B(I,J) = ALPHA*B(I,J)
-  220                     CONTINUE
-                      END IF
+                      DO 220 I = 1,M
+                          B(I,J) = ALPHA*B(I,J)
+  220                 CONTINUE
                      DO 240 K = J + 1,N
-                          IF (A(K,J).NE.ZERO) THEN
-                              DO 230 I = 1,M
-                                  B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  230                         CONTINUE
-                          END IF
+                          DO 230 I = 1,M
+                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
+  230                     CONTINUE
  240                 CONTINUE
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(J,J)
                          DO 250 I = 1,M
-                              B(I,J) = TEMP*B(I,J)
+                              B(I,J) = B(I,J)/A(J,J)
  250                     CONTINUE
                      END IF
  260             CONTINUE
@@ -389,46 +372,34 @@
              IF (UPPER) THEN
                  DO 310 K = N,1,-1
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(K,K)
                          DO 270 I = 1,M
-                              B(I,K) = TEMP*B(I,K)
+                              B(I,K) = B(I,K)/A(K,K)
  270                     CONTINUE
                      END IF
                      DO 290 J = 1,K - 1
-                          IF (A(J,K).NE.ZERO) THEN
-                              TEMP = A(J,K)
-                              DO 280 I = 1,M
-                                  B(I,J) = B(I,J) - TEMP*B(I,K)
-  280                         CONTINUE
-                          END IF
+                          DO 280 I = 1,M
+                              B(I,J) = B(I,J) - A(J,K)*B(I,K)
+  280                     CONTINUE
  290                 CONTINUE
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 300 I = 1,M
-                              B(I,K) = ALPHA*B(I,K)
-  300                     CONTINUE
-                      END IF
+                      DO 300 I = 1,M
+                          B(I,K) = ALPHA*B(I,K)
+  300                 CONTINUE
  310             CONTINUE
              ELSE
                  DO 360 K = 1,N
                      IF (NOUNIT) THEN
-                          TEMP = ONE/A(K,K)
                          DO 320 I = 1,M
-                              B(I,K) = TEMP*B(I,K)
+                              B(I,K) = B(I,K)/A(K,K)
  320                     CONTINUE
                      END IF
                      DO 340 J = K + 1,N
-                          IF (A(J,K).NE.ZERO) THEN
-                              TEMP = A(J,K)
-                              DO 330 I = 1,M
-                                  B(I,J) = B(I,J) - TEMP*B(I,K)
-  330                         CONTINUE
-                          END IF
+                          DO 330 I = 1,M
+                              B(I,J) = B(I,J) - A(J,K)*B(I,K)
+  330                     CONTINUE
  340                 CONTINUE
-                      IF (ALPHA.NE.ONE) THEN
-                          DO 350 I = 1,M
-                              B(I,K) = ALPHA*B(I,K)
-  350                     CONTINUE
-                      END IF
+                      DO 350 I = 1,M
+                          B(I,K) = ALPHA*B(I,K)
+  350                 CONTINUE
  360             CONTINUE
              END IF
          END IF
@@ -140,6 +140,7 @@
 *
 *  =====================================================================
      SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level2 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -154,10 +155,6 @@
 *     ..
 *
 *  =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
      DOUBLE PRECISION TEMP
@@ -224,52 +221,44 @@
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 20 J = N,1,-1
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/A(J,J)
-                          TEMP = X(J)
-                          DO 10 I = J - 1,1,-1
-                              X(I) = X(I) - TEMP*A(I,J)
-   10                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/A(J,J)
+                      TEMP = X(J)
+                      DO 10 I = J - 1,1,-1
+                          X(I) = X(I) - TEMP*A(I,J)
+   10                 CONTINUE
   20             CONTINUE
              ELSE
                  JX = KX + (N-1)*INCX
                  DO 40 J = N,1,-1
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/A(J,J)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 30 I = J - 1,1,-1
-                              IX = IX - INCX
-                              X(IX) = X(IX) - TEMP*A(I,J)
-   30                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 30 I = J - 1,1,-1
+                          IX = IX - INCX
+                          X(IX) = X(IX) - TEMP*A(I,J)
+   30                 CONTINUE
                      JX = JX - INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = 1,N
-                      IF (X(J).NE.ZERO) THEN
-                          IF (NOUNIT) X(J) = X(J)/A(J,J)
-                          TEMP = X(J)
-                          DO 50 I = J + 1,N
-                              X(I) = X(I) - TEMP*A(I,J)
-   50                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(J) = X(J)/A(J,J)
+                      TEMP = X(J)
+                      DO 50 I = J + 1,N
+                          X(I) = X(I) - TEMP*A(I,J)
+   50                 CONTINUE
   60             CONTINUE
              ELSE
                  JX = KX
                  DO 80 J = 1,N
-                      IF (X(JX).NE.ZERO) THEN
-                          IF (NOUNIT) X(JX) = X(JX)/A(J,J)
-                          TEMP = X(JX)
-                          IX = JX
-                          DO 70 I = J + 1,N
-                              IX = IX + INCX
-                              X(IX) = X(IX) - TEMP*A(I,J)
-   70                     CONTINUE
-                      END IF
+                      IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+                      TEMP = X(JX)
+                      IX = JX
+                      DO 70 I = J + 1,N
+                          IX = IX + INCX
+                          X(IX) = X(IX) - TEMP*A(I,J)
+   70                 CONTINUE
                      JX = JX + INCX
   80             CONTINUE
              END IF
@@ -69,6 +69,7 @@
 *>
 *  =====================================================================
      DOUBLE PRECISION FUNCTION DZASUM(N,ZX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -86,11 +86,12 @@
 !> \endverbatim
 !>
 !  =====================================================================
-function DZNRM2( n, x, incx ) 
+function DZNRM2( n, x, incx )
+   implicit none
   integer, parameter :: wp = kind(1.d0)
   real(wp) :: DZNRM2
 !
-!  -- Reference BLAS level1 routine (version 3.9.1) --
+!  -- Reference BLAS level1 routine --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021
@@ -0,0 +1,193 @@
+!> \brief \b ICAMAX
+!
+!  =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+!            http://www.netlib.org/lapack/explore-html/
+!
+!  Definition:
+!  ===========
+!
+!       INTEGER FUNCTION ICAMAX(N,X,INCX)
+!
+!       .. Scalar Arguments ..
+!       INTEGER INCX,N
+!       ..
+!       .. Array Arguments ..
+!       COMPLEX X(*)
+!       ..
+!
+!
+!> \par Purpose:
+!  =============
+!>
+!> \verbatim
+!>
+!>  ICAMAX finds the index of the first element having maximum |Re(.)| + |Im(.)|
+!> \endverbatim
+!
+!  Arguments:
+!  ==========
+!
+!> \param[in] N
+!> \verbatim
+!>          N is INTEGER
+!>         number of elements in input vector(s)
+!> \endverbatim
+!>
+!> \param[in] X
+!> \verbatim
+!>          X is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+!> \endverbatim
+!>
+!> \param[in] INCX
+!> \verbatim
+!>          INCX is INTEGER
+!>         storage spacing between elements of X
+!> \endverbatim
+!
+!  Authors:
+!  ========
+!
+!> James Demmel, University of California Berkeley, USA
+!> Weslley Pereira, National Renewable Energy Laboratory, USA
+!
+!> \ingroup iamax
+!
+!> \par Further Details:
+!  =====================
+!>
+!> \verbatim
+!>
+!> James Demmel et al. Proposed Consistent Exception Handling for the BLAS and
+!> LAPACK, 2022 (https://arxiv.org/abs/2207.09281).
+!>
+!> \endverbatim
+!>
+!  =====================================================================
+integer function icamax(n, x, incx)
+   implicit none
+   integer, parameter :: wp = kind(1.e0)
+!
+!  -- Reference BLAS level1 routine --
+!  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+!
+!  .. Constants ..
+   real(wp), parameter :: hugeval = huge(0.0_wp)
+!
+!  .. Scalar Arguments ..
+   integer :: n, incx
+!
+!  .. Array Arguments ..
+   complex(wp) :: x(*)
+!  ..
+!  .. Local Scalars ..
+   integer :: i, j, ix, jx
+   real(wp) :: val, smax
+   logical :: scaledsmax
+!  ..
+!  .. Intrinsic Functions ..
+   intrinsic :: abs, aimag, huge, real
+!
+!  Quick return if possible
+!
+   icamax = 0
+   if (n < 1 .or. incx < 1) return
+!
+   icamax = 1
+   if (n == 1) return
+!
+   icamax = 0
+   scaledsmax = .false.
+   smax = -1
+!
+!  scaledsmax = .true. indicates that x(icamax) is finite but
+!  abs(real(x(icamax))) + abs(aimag(x(icamax))) overflows
+!
+   if (incx == 1) then
+      ! code for increment equal to 1
+      do i = 1, n
+         if (x(i) /= x(i)) then
+            ! return when first NaN found
+            icamax = i
+            return
+         elseif (abs(real(x(i))) > hugeval .or. abs(aimag(x(i))) > hugeval) then
+            ! keep looking for first NaN
+            do j = i+1, n
+               if (x(j) /= x(j)) then
+                  ! return when first NaN found
+                  icamax = j
+                  return
+               endif
+            enddo
+            ! record location of first Inf
+            icamax = i
+            return
+         else ! still no Inf found yet
+            if (.not. scaledsmax) then
+               ! no abs(real(x(i))) + abs(aimag(x(i))) = Inf yet
+               val = abs(real(x(i))) + abs(aimag(x(i)))
+               if (val > hugeval) then
+                  scaledsmax = .true.
+                  smax = 0.25*abs(real(x(i))) + 0.25*abs(aimag(x(i)))
+                  icamax = i
+               elseif (val > smax) then ! everything finite so far
+                  smax = val
+                  icamax = i
+               endif
+            else ! scaledsmax
+               val = 0.25*abs(real(x(i))) + 0.25*abs(aimag(x(i)))
+               if (val > smax) then
+                  smax = val
+                  icamax = i
+               endif
+            endif
+         endif
+      end do
+   else
+      ! code for increment not equal to 1
+      ix = 1
+      do i = 1, n
+         if (x(ix) /= x(ix)) then
+            ! return when first NaN found
+            icamax = i
+            return
+         elseif (abs(real(x(ix))) > hugeval .or. abs(aimag(x(ix))) > hugeval) then
+            ! keep looking for first NaN
+            jx = ix + incx
+            do j = i+1, n
+               if (x(jx) /= x(jx)) then
+                  ! return when first NaN found
+                  icamax = j
+                  return
+               endif
+               jx = jx + incx
+            enddo
+            ! record location of first Inf
+            icamax = i
+            return
+         else ! still no Inf found yet
+            if (.not. scaledsmax) then
+               ! no abs(real(x(ix))) + abs(aimag(x(ix))) = Inf yet
+               val = abs(real(x(ix))) + abs(aimag(x(ix)))
+               if (val > hugeval) then
+                  scaledsmax = .true.
+                  smax = 0.25*abs(real(x(ix))) + 0.25*abs(aimag(x(ix)))
+                  icamax = i
+               elseif (val > smax) then ! everything finite so far
+                  smax = val
+                  icamax = i
+               endif
+            else ! scaledsmax
+               val = 0.25*abs(real(x(ix))) + 0.25*abs(aimag(x(ix)))
+               if (val > smax) then
+                  smax = val
+                  icamax = i
+               endif
+            endif
+         endif
+         ix = ix + incx
+      end do
+   endif
+end
@@ -68,6 +68,7 @@
 *>
 *  =====================================================================
      INTEGER FUNCTION IDAMAX(N,DX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -68,6 +68,7 @@
 *>
 *  =====================================================================
      INTEGER FUNCTION ISAMAX(N,SX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -0,0 +1,193 @@
+!> \brief \b IZAMAX
+!
+!  =========== DOCUMENTATION ===========
+!
+! Online html documentation available at
+!            http://www.netlib.org/lapack/explore-html/
+!
+!  Definition:
+!  ===========
+!
+!       INTEGER FUNCTION IZAMAX(N,X,INCX)
+!
+!       .. Scalar Arguments ..
+!       INTEGER INCX,N
+!       ..
+!       .. Array Arguments ..
+!       DOUBLE COMPLEX X(*)
+!       ..
+!
+!
+!> \par Purpose:
+!  =============
+!>
+!> \verbatim
+!>
+!>  IZAMAX finds the index of the first element having maximum |Re(.)| + |Im(.)|
+!> \endverbatim
+!
+!  Arguments:
+!  ==========
+!
+!> \param[in] N
+!> \verbatim
+!>          N is INTEGER
+!>         number of elements in input vector(s)
+!> \endverbatim
+!>
+!> \param[in] X
+!> \verbatim
+!>          X is DOUBLE COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+!> \endverbatim
+!>
+!> \param[in] INCX
+!> \verbatim
+!>          INCX is INTEGER
+!>         storage spacing between elements of X
+!> \endverbatim
+!
+!  Authors:
+!  ========
+!
+!> James Demmel, University of California Berkeley, USA
+!> Weslley Pereira, National Renewable Energy Laboratory, USA
+!
+!> \ingroup iamax
+!
+!> \par Further Details:
+!  =====================
+!>
+!> \verbatim
+!>
+!> James Demmel et al. Proposed Consistent Exception Handling for the BLAS and
+!> LAPACK, 2022 (https://arxiv.org/abs/2207.09281).
+!>
+!> \endverbatim
+!>
+!  =====================================================================
+integer function izamax(n, x, incx)
+   implicit none
+   integer, parameter :: wp = kind(1.d0)
+!
+!  -- Reference BLAS level1 routine --
+!  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+!
+!  .. Constants ..
+   real(wp), parameter :: hugeval = huge(0.0_wp)
+!
+!  .. Scalar Arguments ..
+   integer :: n, incx
+!
+!  .. Array Arguments ..
+   complex(wp) :: x(*)
+!  ..
+!  .. Local Scalars ..
+   integer :: i, j, ix, jx
+   real(wp) :: val, smax
+   logical :: scaledsmax
+!  ..
+!  .. Intrinsic Functions ..
+   intrinsic :: abs, dimag, huge, real
+!
+!  Quick return if possible
+!
+   izamax = 0
+   if (n < 1 .or. incx < 1) return
+!
+   izamax = 1
+   if (n == 1) return
+!
+   izamax = 0
+   scaledsmax = .false.
+   smax = -1
+!
+!  scaledsmax = .true. indicates that x(izamax) is finite but
+!  abs(real(x(izamax))) + abs(dimag(x(izamax))) overflows
+!
+   if (incx == 1) then
+      ! code for increment equal to 1
+      do i = 1, n
+         if (x(i) /= x(i)) then
+            ! return when first NaN found
+            izamax = i
+            return
+         elseif (abs(real(x(i))) > hugeval .or. abs(dimag(x(i))) > hugeval) then
+            ! keep looking for first NaN
+            do j = i+1, n
+               if (x(j) /= x(j)) then
+                  ! return when first NaN found
+                  izamax = j
+                  return
+               endif
+            enddo
+            ! record location of first Inf
+            izamax = i
+            return
+         else ! still no Inf found yet
+            if (.not. scaledsmax) then
+               ! no abs(real(x(i))) + abs(dimag(x(i))) = Inf yet
+               val = abs(real(x(i))) + abs(dimag(x(i)))
+               if (val > hugeval) then
+                  scaledsmax = .true.
+                  smax = 0.25*abs(real(x(i))) + 0.25*abs(dimag(x(i)))
+                  izamax = i
+               elseif (val > smax) then ! everything finite so far
+                  smax = val
+                  izamax = i
+               endif
+            else ! scaledsmax
+               val = 0.25*abs(real(x(i))) + 0.25*abs(dimag(x(i)))
+               if (val > smax) then
+                  smax = val
+                  izamax = i
+               endif
+            endif
+         endif
+      end do
+   else
+      ! code for increment not equal to 1
+      ix = 1
+      do i = 1, n
+         if (x(ix) /= x(ix)) then
+            ! return when first NaN found
+            izamax = i
+            return
+         elseif (abs(real(x(ix))) > hugeval .or. abs(dimag(x(ix))) > hugeval) then
+            ! keep looking for first NaN
+            jx = ix + incx
+            do j = i+1, n
+               if (x(jx) /= x(jx)) then
+                  ! return when first NaN found
+                  izamax = j
+                  return
+               endif
+               jx = jx + incx
+            enddo
+            ! record location of first Inf
+            izamax = i
+            return
+         else ! still no Inf found yet
+            if (.not. scaledsmax) then
+               ! no abs(real(x(ix))) + abs(dimag(x(ix))) = Inf yet
+               val = abs(real(x(ix))) + abs(dimag(x(ix)))
+               if (val > hugeval) then
+                  scaledsmax = .true.
+                  smax = 0.25*abs(real(x(ix))) + 0.25*abs(dimag(x(ix)))
+                  izamax = i
+               elseif (val > smax) then ! everything finite so far
+                  smax = val
+                  izamax = i
+               endif
+            else ! scaledsmax
+               val = 0.25*abs(real(x(ix))) + 0.25*abs(dimag(x(ix)))
+               if (val > smax) then
+                  smax = val
+                  izamax = i
+               endif
+            endif
+         endif
+         ix = ix + incx
+      end do
+   endif
+end
@@ -50,6 +50,7 @@
 *
 *  =====================================================================
      LOGICAL FUNCTION LSAME(CA,CB)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -69,6 +69,7 @@
 *>
 *  =====================================================================
      REAL FUNCTION SASUM(N,SX,INCX)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -0,0 +1,148 @@
+*> \brief \b SAXPBY
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SAXPBY(N,SA,SX,INCX,SB,SY,INCY)
+*
+*       .. Scalar Arguments ..
+*       REAL SA,SB
+*       INTEGER INCX,INCY,N
+*       ..
+*       .. Array Arguments ..
+*       REAL SX(*),SY(*)
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*>    SAXPBY constant times a vector plus constant times a vector.
+*>
+*>    Y = ALPHA * X + BETA * Y
+*>
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          number of elements in input vector(s)
+*> \endverbatim
+*>
+*> \param[in] SA
+*> \verbatim
+*>           SA is REAL
+*>           On entry, SA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] SX
+*> \verbatim
+*>          SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          storage spacing between elements of SX
+*> \endverbatim
+*>
+*> \param[in] SB
+*> \verbatim
+*>           SB is REAL
+*>           On entry, SB specifies the scalar beta.
+*> \endverbatim
+*>
+*> \param[in,out] SY
+*> \verbatim
+*>          SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
+*> \endverbatim
+*>
+*> \param[in] INCY
+*> \verbatim
+*>          INCY is INTEGER
+*>         storage spacing between elements of SY
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*> \author Martin Koehler, MPI Magdeburg
+*
+*> \ingroup axpby
+*
+*  =====================================================================
+      SUBROUTINE SAXPBY(N,SA,SX,INCX,SB,SY,INCY)
+      IMPLICIT NONE
+*
+*  -- Reference BLAS level1 routine --
+*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      REAL SA,SB
+      INTEGER INCX,INCY,N
+*     ..
+*     .. Array Arguments ..
+      REAL SX(*),SY(*)
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL SSCAL
+*
+*  =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER I,IX,IY,M,MP1
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC MOD
+*     ..
+      IF (N.LE.0) RETURN
+
+*     Scale if SA.EQ.0
+      IF (SA.EQ.0.0E0 .AND. SB.NE.0.0E0) THEN
+          CALL SSCAL(N, SB, SY, INCY)
+          RETURN
+      END IF
+
+
+      IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
+*
+*        code for both increments equal to 1
+*
+         DO I = 1,N
+            SY(I) = SB*SY(I) + SA*SX(I)
+         END DO
+      ELSE
+*
+*        code for unequal increments or equal increments
+*          not equal to 1
+*
+         IX = 1
+         IY = 1
+         IF (INCX.LT.0) IX = (-N+1)*INCX + 1
+         IF (INCY.LT.0) IY = (-N+1)*INCY + 1
+         DO I = 1,N
+          SY(IY) = SB*SY(IY) + SA*SX(IX)
+          IX = IX + INCX
+          IY = IY + INCY
+         END DO
+      END IF
+      RETURN
+*
+*     End of SAXPBY
+*
+      END
@@ -86,6 +86,7 @@
 *>
 *  =====================================================================
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -43,6 +43,7 @@
 *
 *  =====================================================================
      REAL FUNCTION SCABS1(Z)
+      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
--- a/Show More
+++ b/Show More