add steqr wrappers for comparison

.appveyor.yml

@@ -15,28 +15,22 @@ 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 flang-rt_win-64 cmake ninja
-  - call "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" amd64
+  - 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
   - 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 "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 ..
+  - cmake -G "NMake Makefiles JOM" -DCMAKE_Fortran_COMPILER=flang -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=ON ..
 
 build_script:
   - cmake --build .
 
 test_script:
-  - ctest -j2 --output-on-failure 
+  - ctest -j2

.github/julia/build_tarballs.jl

@@ -1,65 +0,0 @@
-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")

.github/julia/generate_binaries.jl

@@ -1,90 +0,0 @@
-# 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

.github/workflows/cmake.yml

@@ -7,7 +7,6 @@ on:
     - try-github-actions-for-windows
     paths:
     - .github/workflows/cmake.yml
-    - lapack_testing.py
     - '**CMakeLists.txt'
     - 'BLAS/**'
     - 'CBLAS/**'
@@ -22,7 +21,6 @@ on:
   pull_request:
     paths:
     - .github/workflows/cmake.yml
-    - lapack_testing.py
     - '**CMakeLists.txt'
     - 'BLAS/**'
     - 'CBLAS/**'
@@ -50,7 +48,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.
@@ -64,22 +62,25 @@ jobs:
     strategy:
       fail-fast: true
       matrix:
-        os: [ macos-latest, ubuntu-latest, ubuntu-24.04-arm, windows-latest ]
+        os: [ macos-latest, ubuntu-latest, 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: Use GCC-14 on MacOS
+    - name: Install ninja-build tool
+      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3
+
+    - name: Use GCC-11 on MacOS
       if: ${{ matrix.os == 'macos-latest' }}
       run: >
         cmake -B build -G Ninja
-        -D CMAKE_C_COMPILER="gcc-14"
-        -D CMAKE_Fortran_COMPILER="gfortran-14"
+        -D CMAKE_C_COMPILER="gcc-11"
+        -D CMAKE_Fortran_COMPILER="gfortran-11"
         -D USE_FLAT_NAMESPACE:BOOL=ON
 
     - name: Special flags for Windows
@@ -89,8 +90,7 @@ 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 or Ninja.
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
       # 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
+      run: cmake --build build --config ${{env.BUILD_TYPE}}
 
     - name: Test with OpenMP
       working-directory: ${{github.workspace}}/build
@@ -117,127 +117,24 @@ 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
 
-  test-extended-api-only:
-    runs-on: ubuntu-latest
-
-    env:
-      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
-
-    - name: Install ninja-build tool
-      uses: seanmiddleditch/gha-setup-ninja@16b940825621068d98711680b6c3ff92201f8fc0 # v3
-
-    - name: Configure CMake
-      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
-
+    runs-on: ubuntu-latest 
     env:
       BUILD_TYPE: Coverage
       FFLAGS: "-fopenmp"
-
-    steps:
-
+    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.
+      # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
       # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
       run: >
         cmake -B build -G Ninja
@@ -250,75 +147,17 @@ 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: |
-        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 != '' }}
+        echo "Coverage"
+        cmake --build build --target coverage
+        bash <(curl -s https://codecov.io/bash) -X gcov
 
   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
 
@@ -332,12 +171,9 @@ 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=${{env.BUILD_TYPE}}
+        -D CMAKE_BUILD_TYPE=Release
         -D CMAKE_INSTALL_PREFIX=${{github.workspace}}/lapack_install
         -D CBLAS:BOOL=ON
         -D LAPACKE:BOOL=ON
@@ -348,14 +184,10 @@ 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
 
@@ -366,16 +198,13 @@ 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}}
@@ -387,14 +216,12 @@ jobs:
         -D LAPACK_TESTING_USE_PYTHON:BOOL=OFF
 
     - name: Build
-      # 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
+      run: cmake --build build --config ${{env.BUILD_TYPE}}
 
     - name: Test
       working-directory: ${{github.workspace}}/build
       run: |
-        ctest -C ${{env.BUILD_TYPE}} --output-on-failure --schedule-random -j2 -T memcheck > memcheck.out
+        ctest -C ${{env.BUILD_TYPE}} --schedule-random -j2 -T memcheck > memcheck.out
         cat memcheck.out
         if tail -n 1 memcheck.out | grep -q "Memory checking results:"; then
           exit 0
@@ -410,4 +237,4 @@ jobs:
             fi
           done
           exit 0
-        fi
+        fi

.github/workflows/makefile.yml

@@ -90,8 +90,8 @@ jobs:
         echo "DOCSDIR = ${{github.workspace}}/DOCS" >> make.inc
     - name: Alias for GCC compilers
       run: |
-        sudo ln -s $(which gcc-14) /usr/local/bin/gcc
-        sudo ln -s $(which gfortran-14) /usr/local/bin/gfortran
+        sudo ln -s $(which gcc-11) /usr/local/bin/gcc
+        sudo ln -s $(which gfortran-11) /usr/local/bin/gfortran
     - name: Install
       run: |
         make -s -j2 all

.github/workflows/release.yml

@@ -1,239 +0,0 @@
-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 }}

.github/workflows/scorecard.yml

@@ -32,12 +32,12 @@ jobs:
 
     steps:
       - name: "Checkout code"
-        uses: actions/checkout@d632683dd7b4114ad314bca15554477dd762a938 # tag=v4.2.0
+        uses: actions/checkout@c85c95e3d7251135ab7dc9ce3241c5835cc595a9 # v3.5.3
         with:
           persist-credentials: false
 
       - name: "Run analysis"
-        uses: ossf/scorecard-action@62b2cac7ed8198b15735ed49ab1e5cf35480ba46 # v2.4.0
+        uses: ossf/scorecard-action@08b4669551908b1024bb425080c797723083c031 # v2.2.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@b4b15b8c7c6ac21ea08fcf65892d2ee8f75cf882 # v4.4.3
+        uses: actions/upload-artifact@0b7f8abb1508181956e8e162db84b466c27e18ce # v3.1.2
         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@662472033e021d55d94146f66f6058822b0b39fd # v3.27.0
+        uses: github/codeql-action/upload-sarif@f9a7c6738f28efb36e31d49c53a201a9c5d6a476 # v2.14.2
         with:
           sarif_file: results.sarif

.gitignore

@@ -1,9 +1,5 @@
 # ignore objects and archives, anywhere in the tree.
 *.[oa]
-*.so
-*.dll
-*.dylib
-*.pdb
 
 # test in INSTALL
 INSTALL/test*
@@ -27,7 +23,6 @@ CBLAS/examples/cblas_ex1
 CBLAS/examples/cblas_ex2
 
 # LAPACK testing
-/lapack_testing_junit.xml
 TESTING/LIN/xlintst*
 TESTING/EIG/xeigtst*
 TESTING/EIG/xdmd*
@@ -49,6 +44,5 @@ DOCS/man
 DOCS/explore-html
 output_err
 
-# Mod files from compilation in SRC
-SRC/la_constants.mod
-SRC/la_xisnan.mod
+# Editor config files
+.vscode/

BLAS/SRC/CMakeLists.txt

@@ -29,33 +29,31 @@
 #  Level 1 BLAS
 #---------------------------------------------------------
 
-set(LAPACK_INSTALL_EXPORT_NAME ${BLASLIB}-targets)
+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(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(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(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(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(DB1AUX sscal.f isamax.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(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(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.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)
+    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)
 
 #---------------------------------------------------------------------
 #  Auxiliary routines needed by both the Level 2 and Level 3 BLAS
@@ -65,40 +63,34 @@ 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 sskewsymv.f sskewsyr2.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)
 
-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 dskewsymv.f dskewsyr2.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(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 sskewsymm.f
-  sskewsyr2k.f)
+set(SBLAS3 sgemm.f ssymm.f ssyrk.f ssyr2k.f strmm.f strsm.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)
 
-set(DBLAS3
-  dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.f dgemmtr.f dskewsymm.f
-  dskewsyr2k.f)
+set(DBLAS3 dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.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)
 
 
 set(SOURCES)
@@ -117,43 +109,53 @@ if(BUILD_COMPLEX16)
 endif()
 list(REMOVE_DUPLICATES SOURCES)
 
-if(BUILD_DEFAULT_API)
-  add_library(${BLASLIB}_obj OBJECT ${SOURCES})
-endif()
+add_library(${BLASLIB}_obj OBJECT ${SOURCES})
+set_target_properties(${BLASLIB}_obj PROPERTIES POSITION_INDEPENDENT_CODE ON)
 
 if(BUILD_INDEX64_EXT_API)
-  include(ExtendedAPIHelpers)
-  generate_64bit_suffixed_sources(${BLASLIB} SOURCES SOURCES_64)
-
-  add_library(${BLASLIB}_64_obj OBJECT ${SOURCES_64})
+  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})
   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}
-  $<$<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()
+        $<TARGET_OBJECTS:${BLASLIB}_obj>
+        $<$<BOOL:${BUILD_INDEX64_EXT_API}>: $<TARGET_OBJECTS:${BLASLIB}_64_obj>>)
 
 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()

BLAS/SRC/Makefile

@@ -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 saxpby.o scopy.o sdot.o snrm2.o \
+SBLAS1 = isamax.o sasum.o saxpy.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 caxpby.o ccopy.o \
+CBLAS1 = scabs1.o scasum.o scnrm2.o icamax.o caxpy.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 daxpby.o dcopy.o ddot.o dnrm2.o \
+DBLAS1 = idamax.o dasum.o daxpy.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 zaxpby.o zcopy.o \
+ZBLAS1 = dcabs1.o dzasum.o dznrm2.o izamax.o zaxpy.o zcopy.o \
 	zdotc.o zdotu.o zdscal.o zrotg.o zscal.o zswap.o zdrot.o
 $(ZBLAS1): $(FRC)
 
@@ -105,8 +105,7 @@ $(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 \
-	sskewsymv.o sskewsyr2.o
+	sger.o ssyr.o sspr.o ssyr2.o sspr2.o
 $(SBLAS2): $(FRC)
 
 CBLAS2 = cgemv.o cgbmv.o chemv.o chbmv.o chpmv.o \
@@ -116,8 +115,7 @@ $(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 \
-	dskewsymv.o dskewsyr2.o
+	dger.o dsyr.o dspr.o dsyr2.o dspr2.o
 $(DBLAS2): $(FRC)
 
 ZBLAS2 = zgemv.o zgbmv.o zhemv.o zhbmv.o zhpmv.o \
@@ -129,20 +127,18 @@ $(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 \
-	sskewsymm.o sskewsyr2k.o
+SBLAS3 = sgemm.o ssymm.o ssyrk.o ssyr2k.o strmm.o strsm.o
 $(SBLAS3): $(FRC)
 
 CBLAS3 = cgemm.o csymm.o csyrk.o csyr2k.o ctrmm.o ctrsm.o \
-	chemm.o cherk.o cher2k.o cgemmtr.o
+	chemm.o cherk.o cher2k.o
 $(CBLAS3): $(FRC)
 
-DBLAS3 = dgemm.o dsymm.o dsyrk.o dsyr2k.o dtrmm.o dtrsm.o dgemmtr.o \
-	dskewsymm.o dskewsyr2k.o
+DBLAS3 = dgemm.o dsymm.o dsyrk.o dsyr2k.o dtrmm.o dtrsm.o
 $(DBLAS3): $(FRC)
 
 ZBLAS3 = zgemm.o zsymm.o zsyrk.o zsyr2k.o ztrmm.o ztrsm.o \
-	zhemm.o zherk.o zher2k.o zgemmtr.o
+	zhemm.o zherk.o zher2k.o
 $(ZBLAS3): $(FRC)
 
 ALLOBJ = $(SBLAS1) $(SBLAS2) $(SBLAS3) $(DBLAS1) $(DBLAS2) $(DBLAS3) \

BLAS/SRC/caxpby.f

@@ -1,144 +0,0 @@
-*> \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

BLAS/SRC/caxpy.f

@@ -85,7 +85,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -103,16 +102,13 @@
 *
 *     .. Local Scalars ..
       INTEGER I,IX,IY
-      COMPLEX CDUM
 *     ..
-*     .. Statement Functions ..
-      REAL CABS1
-*     ..
-*     .. Statement Function definitions ..
-      CABS1(CDUM) = ABS(REAL(CDUM)) + ABS(AIMAG(CDUM))
+*     .. External Functions ..
+      REAL SCABS1
+      EXTERNAL SCABS1
 *     ..
       IF (N.LE.0) RETURN
-      IF (CABS1(CA).EQ.0.0E+0) RETURN
+      IF (SCABS1(CA).EQ.0.0E+0) RETURN
       IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
 *
 *        code for both increments equal to 1

BLAS/SRC/ccopy.f

@@ -78,7 +78,6 @@
 *>
 *  =====================================================================
       SUBROUTINE CCOPY(N,CX,INCX,CY,INCY)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/cdotc.f

@@ -80,7 +80,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cdotu.f

@@ -80,7 +80,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cgbmv.f

@@ -187,7 +187,6 @@
 *  =====================================================================
       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,    --

BLAS/SRC/cgemm.f

@@ -35,16 +35,6 @@
 *>
 *> 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:
@@ -102,9 +92,7 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is COMPLEX
-*>           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.
+*>           On entry, ALPHA specifies the scalar alpha.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -114,10 +102,7 @@
 *>           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, 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.
+*>           matrix A.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -136,10 +121,7 @@
 *>           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, 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.
+*>           matrix B.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -154,19 +136,16 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is COMPLEX
-*>           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.
+*>           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 COMPLEX array, dimension ( LDC, N )
 *>           Before entry, the leading  m by n  part of the array  C must
-*>           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.
+*>           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  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -206,7 +185,6 @@
 *  =====================================================================
       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,    --

BLAS/SRC/cgemmtr.f

@@ -1,569 +0,0 @@
-*> \brief \b CGEMMTR
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at
-*            http://www.netlib.org/lapack/explore-html/
-*
-*  Definition:
-*  ===========
-*
-*       SUBROUTINE CGEMMTR(UPLO,TRANSA,TRANSB,N,K,ALPHA,A,LDA,B,LDB,BETA,
-*                         C,LDC)
-*
-*       .. Scalar Arguments ..
-*       COMPLEX ALPHA,BETA
-*       INTEGER K,LDA,LDB,LDC,N
-*       CHARACTER TRANSA,TRANSB, UPLO
-*       ..
-*       .. Array Arguments ..
-*       COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
-*       ..
-*
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> CGEMMTR  performs one of the matrix-matrix operations
-*>
-*>    C := alpha*op( A )*op( B ) + beta*C,
-*>
-*> where  op( X ) is one of
-*>
-*>    op( X ) = X   or   op( X ) = X**T,
-*>
-*> alpha and beta are scalars, and A, B and C are matrices, with op( A )
-*> an n by k matrix,  op( B )  a  k by n matrix and  C an n by n matrix.
-*> Thereby, the routine only accesses and updates the upper or lower
-*> triangular part of the result matrix C. This behaviour can be used if
-*> the resulting matrix C is known to be Hermitian or symmetric.
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] UPLO
-*> \verbatim
-*>          UPLO is CHARACTER*1
-*>           On entry, UPLO specifies whether the lower or the upper
-*>           triangular part of C is access and updated.
-*>
-*>              UPLO = 'L' or 'l', the lower triangular part of C is used.
-*>
-*>              UPLO = 'U' or 'u', the upper triangular part of C is used.
-*> \endverbatim
-*
-*> \param[in] TRANSA
-*> \verbatim
-*>          TRANSA is CHARACTER*1
-*>           On entry, TRANSA specifies the form of op( A ) to be used in
-*>           the matrix multiplication as follows:
-*>
-*>              TRANSA = 'N' or 'n',  op( A ) = A.
-*>
-*>              TRANSA = 'T' or 't',  op( A ) = A**T.
-*>
-*>              TRANSA = 'C' or 'c',  op( A ) = A**H.
-*> \endverbatim
-*>
-*> \param[in] TRANSB
-*> \verbatim
-*>          TRANSB is CHARACTER*1
-*>           On entry, TRANSB specifies the form of op( B ) to be used in
-*>           the matrix multiplication as follows:
-*>
-*>              TRANSB = 'N' or 'n',  op( B ) = B.
-*>
-*>              TRANSB = 'T' or 't',  op( B ) = B**T.
-*>
-*>              TRANSB = 'C' or 'c',  op( B ) = B**H.
-*> \endverbatim
-*>
-*> \param[in] N
-*> \verbatim
-*>          N is INTEGER
-*>           On entry,  N specifies the number of rows and columns of
-*>           the matrix C, the number of columns of op(B) and the number
-*>           of rows of op(A).  N must be at least zero.
-*> \endverbatim
-*>
-*> \param[in] K
-*> \verbatim
-*>          K is INTEGER
-*>           On entry,  K  specifies  the number of columns of the matrix
-*>           op( A ) and the number of rows of the matrix op( B ). K must
-*>           be at least  zero.
-*> \endverbatim
-*>
-*> \param[in] ALPHA
-*> \verbatim
-*>          ALPHA is COMPLEX.
-*>           On entry, ALPHA specifies the scalar alpha.
-*> \endverbatim
-*>
-*> \param[in] A
-*> \verbatim
-*>          A is COMPLEX array, dimension ( LDA, ka ), where ka is
-*>           k  when  TRANSA = 'N' or 'n',  and is  n  otherwise.
-*>           Before entry with  TRANSA = 'N' or 'n',  the leading  n 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.
-*> \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  TRANSA = '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 COMPLEX array, dimension ( LDB, kb ), where kb is
-*>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
-*>           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.
-*> \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  TRANSB = 'N' or 'n' then
-*>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
-*>           least  max( 1, n ).
-*> \endverbatim
-*>
-*> \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.
-*> \endverbatim
-*>
-*> \param[in,out] C
-*> \verbatim
-*>          C is COMPLEX array, dimension ( LDC, N )
-*>           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 triangular part of the matrix
-*>           C  is overwritten by the n by n matrix
-*>           ( alpha*op( A )*op( B ) + beta*C ).
-*> \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 Martin Koehler
-*
-*> \ingroup gemmtr
-*
-*> \par Further Details:
-*  =====================
-*>
-*> \verbatim
-*>
-*>  Level 3 Blas routine.
-*>
-*>  -- Written on 19-July-2023.
-*>     Martin Koehler, MPI Magdeburg
-*> \endverbatim
-*>
-*  =====================================================================
-      SUBROUTINE CGEMMTR(UPLO,TRANSA,TRANSB,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 ..
-      COMPLEX ALPHA,BETA
-      INTEGER K,LDA,LDB,LDC,N
-      CHARACTER TRANSA,TRANSB,UPLO
-*     ..
-*     .. Array Arguments ..
-      COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
-*     ..
-*
-*  =====================================================================
-*
-*     .. External Functions ..
-      LOGICAL LSAME
-      EXTERNAL LSAME
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL XERBLA
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC CONJG,MAX
-*     ..
-*     .. Local Scalars ..
-      COMPLEX TEMP
-      INTEGER I,INFO,J,L,NROWA,NROWB,ISTART, ISTOP
-      LOGICAL CONJA,CONJB,NOTA,NOTB,UPPER
-*     ..
-*     .. Parameters ..
-      COMPLEX ONE
-      PARAMETER (ONE= (1.0E+0,0.0E+0))
-      COMPLEX ZERO
-      PARAMETER (ZERO= (0.0E+0,0.0E+0))
-*     ..
-*
-*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
-*     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
-*     B  respectively are to be  transposed but  not conjugated  and set
-*     NROWA and  NROWB  as the number of rows of  A  and  B  respectively.
-*
-      NOTA = LSAME(TRANSA,'N')
-      NOTB = LSAME(TRANSB,'N')
-      CONJA = LSAME(TRANSA,'C')
-      CONJB = LSAME(TRANSB,'C')
-      IF (NOTA) THEN
-          NROWA = N
-      ELSE
-          NROWA = K
-      END IF
-      IF (NOTB) THEN
-          NROWB = K
-      ELSE
-          NROWB = N
-      END IF
-      UPPER = LSAME(UPLO, 'U')
-
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF ((.NOT. UPPER) .AND. (.NOT. LSAME(UPLO, 'L'))) THEN
-          INFO = 1
-      ELSE IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND.
-     +    (.NOT.LSAME(TRANSA,'T'))) THEN
-          INFO = 2
-      ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND.
-     +         (.NOT.LSAME(TRANSB,'T'))) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (K.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
-          INFO = 8
-      ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
-          INFO = 10
-      ELSE IF (LDC.LT.MAX(1,N)) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('CGEMMTR',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-*     And when  alpha.eq.zero.
-*
-      IF (ALPHA.EQ.ZERO) THEN
-          IF (BETA.EQ.ZERO) THEN
-              DO 20 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 10 I = ISTART, ISTOP
-                      C(I,J) = ZERO
-   10             CONTINUE
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-                  DO 30 I = ISTART, ISTOP
-                      C(I,J) = BETA*C(I,J)
-   30             CONTINUE
-   40         CONTINUE
-          END IF
-          RETURN
-      END IF
-*
-*     Start the operations.
-*
-      IF (NOTB) THEN
-          IF (NOTA) THEN
-*
-*           Form  C := alpha*A*B + beta*C.
-*
-              DO 90 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-                  IF (BETA.EQ.ZERO) THEN
-                      DO 50 I = ISTART, ISTOP
-                          C(I,J) = ZERO
-   50                 CONTINUE
-                  ELSE IF (BETA.NE.ONE) THEN
-                      DO 60 I = ISTART, ISTOP
-                          C(I,J) = BETA*C(I,J)
-   60                 CONTINUE
-                  END IF
-                  DO 80 L = 1,K
-                      TEMP = ALPHA*B(L,J)
-                      DO 70 I = ISTART, ISTOP
-                          C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                 CONTINUE
-   80             CONTINUE
-   90         CONTINUE
-          ELSE IF (CONJA) THEN
-*
-*           Form  C := alpha*A**H*B + beta*C.
-*
-              DO 120 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 110 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 100 L = 1,K
-                          TEMP = TEMP + CONJG(A(L,I))*B(L,J)
-  100                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  110             CONTINUE
-  120         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A**T*B + beta*C
-*
-              DO 150 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 140 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 130 L = 1,K
-                          TEMP = TEMP + A(L,I)*B(L,J)
-  130                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  140             CONTINUE
-  150         CONTINUE
-          END IF
-      ELSE IF (NOTA) THEN
-          IF (CONJB) THEN
-*
-*           Form  C := alpha*A*B**H + beta*C.
-*
-              DO 200 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  IF (BETA.EQ.ZERO) THEN
-                      DO 160 I = ISTART,ISTOP
-                          C(I,J) = ZERO
-  160                 CONTINUE
-                  ELSE IF (BETA.NE.ONE) THEN
-                      DO 170 I = ISTART, ISTOP
-                          C(I,J) = BETA*C(I,J)
-  170                 CONTINUE
-                  END IF
-                  DO 190 L = 1,K
-                      TEMP = ALPHA*CONJG(B(J,L))
-                      DO 180 I = ISTART, ISTOP
-                          C(I,J) = C(I,J) + TEMP*A(I,L)
-  180                 CONTINUE
-  190             CONTINUE
-  200         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A*B**T + beta*C
-*
-              DO 250 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  IF (BETA.EQ.ZERO) THEN
-                      DO 210 I = ISTART, ISTOP
-                          C(I,J) = ZERO
-  210                 CONTINUE
-                  ELSE IF (BETA.NE.ONE) THEN
-                      DO 220 I = ISTART, ISTOP
-                          C(I,J) = BETA*C(I,J)
-  220                 CONTINUE
-                  END IF
-                  DO 240 L = 1,K
-                      TEMP = ALPHA*B(J,L)
-                      DO 230 I = ISTART, ISTOP
-                          C(I,J) = C(I,J) + TEMP*A(I,L)
-  230                 CONTINUE
-  240             CONTINUE
-  250         CONTINUE
-          END IF
-      ELSE IF (CONJA) THEN
-          IF (CONJB) THEN
-*
-*           Form  C := alpha*A**H*B**H + beta*C.
-*
-              DO 280 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 270 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 260 L = 1,K
-                          TEMP = TEMP + CONJG(A(L,I))*CONJG(B(J,L))
-  260                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  270             CONTINUE
-  280         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A**H*B**T + beta*C
-*
-              DO 310 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 300 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 290 L = 1,K
-                          TEMP = TEMP + CONJG(A(L,I))*B(J,L)
-  290                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  300             CONTINUE
-  310         CONTINUE
-          END IF
-      ELSE
-          IF (CONJB) THEN
-*
-*           Form  C := alpha*A**T*B**H + beta*C
-*
-              DO 340 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 330 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 320 L = 1,K
-                          TEMP = TEMP + A(L,I)*CONJG(B(J,L))
-  320                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  330             CONTINUE
-  340         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A**T*B**T + beta*C
-*
-              DO 370 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 360 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 350 L = 1,K
-                          TEMP = TEMP + A(L,I)*B(J,L)
-  350                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  360             CONTINUE
-  370         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of CGEMMTR
-*
-      END

BLAS/SRC/cgemv.f

@@ -157,7 +157,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cgerc.f

@@ -127,7 +127,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cgeru.f

@@ -127,7 +127,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chbmv.f

@@ -184,7 +184,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chemm.f

@@ -188,7 +188,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chemv.f

@@ -151,7 +151,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cher.f

@@ -132,7 +132,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cher2.f

@@ -147,7 +147,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cher2k.f

@@ -194,7 +194,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/cherk.f

@@ -170,7 +170,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chpmv.f

@@ -146,7 +146,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chpr.f

@@ -127,7 +127,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/chpr2.f

@@ -142,7 +142,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/crotg.f90

@@ -86,7 +86,6 @@
 !
 !  =====================================================================
 subroutine CROTG( a, b, c, s )
-   implicit none
    integer, parameter :: wp = kind(1.e0)
 !
 !  -- Reference BLAS level1 routine --

BLAS/SRC/cscal.f

@@ -75,7 +75,6 @@
 *>
 *  =====================================================================
       SUBROUTINE CSCAL(N,CA,CX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/csrot.f

@@ -95,7 +95,6 @@
 *
 *  =====================================================================
       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,    --

BLAS/SRC/csscal.f

@@ -75,7 +75,6 @@
 *>
 *  =====================================================================
       SUBROUTINE CSSCAL(N,SA,CX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/cswap.f

@@ -78,7 +78,6 @@
 *>
 *  =====================================================================
       SUBROUTINE CSWAP(N,CX,INCX,CY,INCY)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/csymm.f

@@ -186,7 +186,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/csyr2k.f

@@ -185,7 +185,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/csyrk.f

@@ -164,7 +164,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/ctbmv.f

@@ -183,7 +183,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -198,6 +197,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -268,24 +271,28 @@
               KPLUS1 = K + 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
                       IF (J.GT.K) KX = KX + INCX
    40             CONTINUE
@@ -293,25 +300,29 @@
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
                       IF ((N-J).GE.K) KX = KX - INCX
    80             CONTINUE

BLAS/SRC/ctbsv.f

@@ -186,7 +186,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -201,6 +200,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -271,51 +274,59 @@
               KPLUS1 = K + 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
    20             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 40 J = N,1,-1
                       KX = KX - INCX
-                      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
+                      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
                       JX = JX - INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
                       KX = KX + INCX
-                      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
+                      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
                       JX = JX + INCX
    80             CONTINUE
               END IF

BLAS/SRC/ctpmv.f

@@ -139,7 +139,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -154,6 +153,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -220,25 +223,29 @@
               KK = 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
                       KK = KK + J
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
                       KK = KK + J
    40             CONTINUE
@@ -247,26 +254,30 @@
               KK = (N* (N+1))/2
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
                       KK = KK - (N-J+1)
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
                       KK = KK - (N-J+1)
    80             CONTINUE

BLAS/SRC/ctpsv.f

@@ -141,7 +141,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -156,6 +155,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -222,25 +225,29 @@
               KK = (N* (N+1))/2
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
                       KK = KK - J
    20             CONTINUE
               ELSE
                   JX = KX + (N-1)*INCX
                   DO 40 J = N,1,-1
-                      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
+                      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
                       JX = JX - INCX
                       KK = KK - J
    40             CONTINUE
@@ -249,25 +256,29 @@
               KK = 1
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
                       KK = KK + (N-J+1)
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
-                      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
+                      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
                       JX = JX + INCX
                       KK = KK + (N-J+1)
    80             CONTINUE

BLAS/SRC/ctrmm.f

@@ -174,7 +174,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -276,23 +275,27 @@
               IF (UPPER) THEN
                   DO 50 J = 1,N
                       DO 40 K = 1,M
-                          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
+                          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
    40                 CONTINUE
    50             CONTINUE
               ELSE
                   DO 80 J = 1,N
                       DO 70 K = M,1,-1
-                          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
+                          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
    70                 CONTINUE
    80             CONTINUE
               END IF
@@ -351,10 +354,12 @@
                           B(I,J) = TEMP*B(I,J)
   170                 CONTINUE
                       DO 190 K = 1,J - 1
-                          TEMP = ALPHA*A(K,J)
-                          DO 180 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  180                     CONTINUE
+                          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
   190                 CONTINUE
   200             CONTINUE
               ELSE
@@ -365,10 +370,12 @@
                           B(I,J) = TEMP*B(I,J)
   210                 CONTINUE
                       DO 230 K = J + 1,N
-                          TEMP = ALPHA*A(K,J)
-                          DO 220 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  220                     CONTINUE
+                          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
   230                 CONTINUE
   240             CONTINUE
               END IF
@@ -379,14 +386,16 @@
               IF (UPPER) THEN
                   DO 280 K = 1,N
                       DO 260 J = 1,K - 1
-                          IF (NOCONJ) THEN
-                              TEMP = ALPHA*A(J,K)
-                          ELSE
-                              TEMP = ALPHA*CONJG(A(J,K))
+                          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
                           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
@@ -405,14 +414,16 @@
               ELSE
                   DO 320 K = N,1,-1
                       DO 300 J = K + 1,N
-                          IF (NOCONJ) THEN
-                              TEMP = ALPHA*A(J,K)
-                          ELSE
-                              TEMP = ALPHA*CONJG(A(J,K))
+                          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
                           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

BLAS/SRC/ctrmv.f

@@ -144,7 +144,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -159,6 +158,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -226,45 +229,53 @@
           IF (LSAME(UPLO,'U')) THEN
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
    80             CONTINUE
               END IF

BLAS/SRC/ctrsm.f

@@ -177,7 +177,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -210,6 +209,8 @@
       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))
 *     ..
@@ -276,26 +277,34 @@
 *
               IF (UPPER) THEN
                   DO 60 J = 1,N
-                       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
+                      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
    50                 CONTINUE
    60             CONTINUE
-               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
+              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
    90                 CONTINUE
   100             CONTINUE
               END IF
@@ -349,33 +358,43 @@
 *
               IF (UPPER) THEN
                   DO 230 J = 1,N
-                      DO 190 I = 1,M
-                          B(I,J) = ALPHA*B(I,J)
-  190                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 190 I = 1,M
+                              B(I,J) = ALPHA*B(I,J)
+  190                     CONTINUE
+                      END IF
                       DO 210 K = 1,J - 1
-                          DO 200 I = 1,M
-                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  200                     CONTINUE
+                          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
   210                 CONTINUE
                       IF (NOUNIT) THEN
+                          TEMP = ONE/A(J,J)
                           DO 220 I = 1,M
-                              B(I,J) = B(I,J)/A(J,J)
+                              B(I,J) = TEMP*B(I,J)
   220                     CONTINUE
                       END IF
   230             CONTINUE
               ELSE
                   DO 280 J = N,1,-1
-                      DO 240 I = 1,M
-                          B(I,J) = ALPHA*B(I,J)
-  240                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 240 I = 1,M
+                              B(I,J) = ALPHA*B(I,J)
+  240                     CONTINUE
+                      END IF
                       DO 260 K = J + 1,N
-                          DO 250 I = 1,M
-                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  250                     CONTINUE
+                          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
   260                 CONTINUE
                       IF (NOUNIT) THEN
+                          TEMP = ONE/A(J,J)
                           DO 270 I = 1,M
-                              B(I,J) = B(I,J)/A(J,J)
+                              B(I,J) = TEMP*B(I,J)
   270                     CONTINUE
                       END IF
   280             CONTINUE
@@ -389,55 +408,61 @@
                   DO 330 K = N,1,-1
                       IF (NOUNIT) THEN
                           IF (NOCONJ) THEN
-                              DO 290 I = 1,M
-                                  B(I,K) = B(I,K)/A(K,K)
-  290                         CONTINUE
+                              TEMP = ONE/A(K,K)
                           ELSE
-                              DO 390 I = 1,M
-                                  B(I,K) = B(I,K)/CONJG(A(K,K))
-  390                         CONTINUE
+                              TEMP = ONE/CONJG(A(K,K))
                           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 (NOCONJ) THEN
-                              TEMP = A(J,K)
-                          ELSE
-                              TEMP = CONJG(A(J,K))
+                          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
                           END IF
-                          DO 300 I = 1,M
-                              B(I,J) = B(I,J) - TEMP*B(I,K)
-  300                     CONTINUE
   310                 CONTINUE
-                      DO 320 I = 1,M
-                          B(I,K) = ALPHA*B(I,K)
-  320                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 320 I = 1,M
+                              B(I,K) = ALPHA*B(I,K)
+  320                     CONTINUE
+                      END IF
   330             CONTINUE
               ELSE
                   DO 380 K = 1,N
                       IF (NOUNIT) THEN
                           IF (NOCONJ) THEN
-                              DO 340 I = 1,M
-                                  B(I,K) = B(I,K)/A(K,K)
-  340                         CONTINUE
+                              TEMP = ONE/A(K,K)
                           ELSE
-                              DO 400 I = 1,M
-                                  B(I,K) = B(I,K)/CONJG(A(K,K))
-  400                         CONTINUE
+                              TEMP = ONE/CONJG(A(K,K))
                           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 (NOCONJ) THEN
-                              TEMP = A(J,K)
-                          ELSE
-                              TEMP = CONJG(A(J,K))
+                          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
                           END IF
-                          DO 350 I = 1,M
-                              B(I,J) = B(I,J) - TEMP*B(I,K)
-  350                     CONTINUE
   360                 CONTINUE
-                      DO 370 I = 1,M
-                          B(I,K) = ALPHA*B(I,K)
-  370                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 370 I = 1,M
+                              B(I,K) = ALPHA*B(I,K)
+  370                     CONTINUE
+                      END IF
   380             CONTINUE
               END IF
           END IF

BLAS/SRC/ctrsv.f

@@ -146,7 +146,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -161,6 +160,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX ZERO
+      PARAMETER (ZERO= (0.0E+0,0.0E+0))
 *     ..
 *     .. Local Scalars ..
       COMPLEX TEMP
@@ -228,44 +231,52 @@
           IF (LSAME(UPLO,'U')) THEN
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
    20             CONTINUE
               ELSE
                   JX = KX + (N-1)*INCX
                   DO 40 J = N,1,-1
-                      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
+                      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
                       JX = JX - INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
-                      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
+                      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
                       JX = JX + INCX
    80             CONTINUE
               END IF

BLAS/SRC/dasum.f

@@ -68,7 +68,6 @@
 *>
 *  =====================================================================
       DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/daxpby.f

@@ -1,149 +0,0 @@
-*> \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

BLAS/SRC/daxpy.f

@@ -86,7 +86,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dcabs1.f

@@ -44,7 +44,6 @@
 *
 *  =====================================================================
       DOUBLE PRECISION FUNCTION DCABS1(Z)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/dcopy.f

@@ -79,7 +79,6 @@
 *>
 *  =====================================================================
       SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/ddot.f

@@ -79,7 +79,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dgbmv.f

@@ -185,7 +185,6 @@
 *  =====================================================================
       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,    --

BLAS/SRC/dgemm.f

@@ -35,16 +35,6 @@
 *>
 *> 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:
@@ -61,9 +51,6 @@
 *>              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
@@ -77,9 +64,6 @@
 *>              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
@@ -108,9 +92,7 @@
 *> \param[in] ALPHA
 *> \verbatim
 *>          ALPHA is DOUBLE PRECISION.
-*>           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.
+*>           On entry, ALPHA specifies the scalar alpha.
 *> \endverbatim
 *>
 *> \param[in] A
@@ -120,10 +102,7 @@
 *>           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, 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.
+*>           matrix A.
 *> \endverbatim
 *>
 *> \param[in] LDA
@@ -142,10 +121,7 @@
 *>           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, 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.
+*>           matrix B.
 *> \endverbatim
 *>
 *> \param[in] LDB
@@ -160,19 +136,16 @@
 *> \param[in] BETA
 *> \verbatim
 *>          BETA is DOUBLE PRECISION.
-*>           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.
+*>           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 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.
+*>           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  matrix
 *>           ( alpha*op( A )*op( B ) + beta*C ).
 *> \endverbatim
@@ -212,7 +185,6 @@
 *  =====================================================================
       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,    --

BLAS/SRC/dgemmtr.f

@@ -1,431 +0,0 @@
-*> \brief \b DGEMMTR
-*
-*  =========== DOCUMENTATION ===========
-*
-* Online html documentation available at
-*            http://www.netlib.org/lapack/explore-html/
-*
-*  Definition:
-*  ===========
-*
-*       SUBROUTINE DGEMMTR(UPLO,TRANSA,TRANSB,N,K,ALPHA,A,LDA,B,LDB,BETA,
-*                         C,LDC)
-*
-*       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA,BETA
-*       INTEGER K,LDA,LDB,LDC,N
-*       CHARACTER TRANSA,TRANSB, UPLO
-*       ..
-*       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
-*       ..
-*
-*
-*> \par Purpose:
-*  =============
-*>
-*> \verbatim
-*>
-*> DGEMMTR  performs one of the matrix-matrix operations
-*>
-*>    C := alpha*op( A )*op( B ) + beta*C,
-*>
-*> where  op( X ) is one of
-*>
-*>    op( X ) = X   or   op( X ) = X**T,
-*>
-*> alpha and beta are scalars, and A, B and C are matrices, with op( A )
-*> an n by k matrix,  op( B )  a  k by n matrix and  C an n by n matrix.
-*> Thereby, the routine only accesses and updates the upper or lower
-*> triangular part of the result matrix C. This behaviour can be used if
-*> the resulting matrix C is known to be symmetric.
-*> \endverbatim
-*
-*  Arguments:
-*  ==========
-*
-*> \param[in] UPLO
-*> \verbatim
-*>          UPLO is CHARACTER*1
-*>           On entry, UPLO specifies whether the lower or the upper
-*>           triangular part of C is access and updated.
-*>
-*>              UPLO = 'L' or 'l', the lower triangular part of C is used.
-*>
-*>              UPLO = 'U' or 'u', the upper triangular part of C is used.
-*> \endverbatim
-*
-*> \param[in] TRANSA
-*> \verbatim
-*>          TRANSA is CHARACTER*1
-*>           On entry, TRANSA specifies the form of op( A ) to be used in
-*>           the matrix multiplication as follows:
-*>
-*>              TRANSA = 'N' or 'n',  op( A ) = A.
-*>
-*>              TRANSA = 'T' or 't',  op( A ) = A**T.
-*>
-*>              TRANSA = 'C' or 'c',  op( A ) = A**T.
-*> \endverbatim
-*>
-*> \param[in] TRANSB
-*> \verbatim
-*>          TRANSB is CHARACTER*1
-*>           On entry, TRANSB specifies the form of op( B ) to be used in
-*>           the matrix multiplication as follows:
-*>
-*>              TRANSB = 'N' or 'n',  op( B ) = B.
-*>
-*>              TRANSB = 'T' or 't',  op( B ) = B**T.
-*>
-*>              TRANSB = 'C' or 'c',  op( B ) = B**T.
-*> \endverbatim
-*>
-*> \param[in] N
-*> \verbatim
-*>          N is INTEGER
-*>           On entry,  N specifies the number of rows and columns of
-*>           the matrix C, the number of columns of op(B) and the number
-*>           of rows of op(A).  N must be at least zero.
-*> \endverbatim
-*>
-*> \param[in] K
-*> \verbatim
-*>          K is INTEGER
-*>           On entry,  K  specifies  the number of columns of the matrix
-*>           op( A ) and the number of rows of the matrix op( 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  TRANSA = 'N' or 'n',  and is  n  otherwise.
-*>           Before entry with  TRANSA = 'N' or 'n',  the leading  n 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.
-*> \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  TRANSA = '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
-*>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
-*>           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.
-*> \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  TRANSB = 'N' or 'n' then
-*>           LDB must be at least  max( 1, k ), otherwise  LDB must be at
-*>           least  max( 1, n ).
-*> \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  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 triangular part of the matrix
-*>           C  is overwritten by the n by n matrix
-*>           ( alpha*op( A )*op( B ) + beta*C ).
-*> \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 Martin Koehler
-*
-*> \ingroup gemmtr
-*
-*> \par Further Details:
-*  =====================
-*>
-*> \verbatim
-*>
-*>  Level 3 Blas routine.
-*>
-*>  -- Written on 19-July-2023.
-*>     Martin Koehler, MPI Magdeburg
-*> \endverbatim
-*>
-*  =====================================================================
-      SUBROUTINE DGEMMTR(UPLO,TRANSA,TRANSB,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 TRANSA,TRANSB,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 TEMP
-      INTEGER I,INFO,J,L,NROWA,NROWB, ISTART, ISTOP
-      LOGICAL NOTA,NOTB, UPPER
-*     ..
-*     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
-      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
-*     ..
-*
-*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
-*     transposed and set  NROWA and NROWB  as the number of rows of  A
-*     and  B  respectively.
-*
-      NOTA = LSAME(TRANSA,'N')
-      NOTB = LSAME(TRANSB,'N')
-      IF (NOTA) THEN
-          NROWA = N
-      ELSE
-          NROWA = K
-      END IF
-      IF (NOTB) THEN
-          NROWB = K
-      ELSE
-          NROWB = N
-      END IF
-      UPPER = LSAME(UPLO, 'U')
-*
-*     Test the input parameters.
-*
-      INFO = 0
-      IF ((.NOT. UPPER) .AND. (.NOT. LSAME(UPLO, 'L'))) THEN
-          INFO = 1
-      ELSE IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
-     +    (.NOT.LSAME(TRANSA,'T'))) THEN
-          INFO = 2
-      ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
-     +         (.NOT.LSAME(TRANSB,'T'))) THEN
-          INFO = 3
-      ELSE IF (N.LT.0) THEN
-          INFO = 4
-      ELSE IF (K.LT.0) THEN
-          INFO = 5
-      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
-          INFO = 8
-      ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
-          INFO = 10
-      ELSE IF (LDC.LT.MAX(1,N)) THEN
-          INFO = 13
-      END IF
-      IF (INFO.NE.0) THEN
-          CALL XERBLA('DGEMMTR',INFO)
-          RETURN
-      END IF
-*
-*     Quick return if possible.
-*
-      IF (N.EQ.0) RETURN
-*
-*     And if  alpha.eq.zero.
-*
-      IF (ALPHA.EQ.ZERO) THEN
-          IF (BETA.EQ.ZERO) THEN
-              DO 20 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 10 I = ISTART, ISTOP
-                      C(I,J) = ZERO
-   10             CONTINUE
-   20         CONTINUE
-          ELSE
-              DO 40 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 30 I = ISTART, ISTOP
-                      C(I,J) = BETA*C(I,J)
-   30             CONTINUE
-   40         CONTINUE
-          END IF
-          RETURN
-      END IF
-*
-*     Start the operations.
-*
-      IF (NOTB) THEN
-          IF (NOTA) THEN
-*
-*           Form  C := alpha*A*B + beta*C.
-*
-              DO 90 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-                  IF (BETA.EQ.ZERO) THEN
-                      DO 50 I = ISTART, ISTOP
-                          C(I,J) = ZERO
-   50                 CONTINUE
-                  ELSE IF (BETA.NE.ONE) THEN
-                      DO 60 I = ISTART, ISTOP
-                          C(I,J) = BETA*C(I,J)
-   60                 CONTINUE
-                  END IF
-                  DO 80 L = 1,K
-                      TEMP = ALPHA*B(L,J)
-                      DO 70 I = ISTART, ISTOP
-                          C(I,J) = C(I,J) + TEMP*A(I,L)
-   70                 CONTINUE
-   80             CONTINUE
-   90         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A**T*B + beta*C
-*
-              DO 120 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 110 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 100 L = 1,K
-                          TEMP = TEMP + A(L,I)*B(L,J)
-  100                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  110             CONTINUE
-  120         CONTINUE
-          END IF
-      ELSE
-          IF (NOTA) THEN
-*
-*           Form  C := alpha*A*B**T + beta*C
-*
-              DO 170 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  IF (BETA.EQ.ZERO) THEN
-                      DO 130 I = ISTART,ISTOP
-                          C(I,J) = ZERO
-  130                 CONTINUE
-                  ELSE IF (BETA.NE.ONE) THEN
-                      DO 140 I = ISTART,ISTOP
-                          C(I,J) = BETA*C(I,J)
-  140                 CONTINUE
-                  END IF
-                  DO 160 L = 1,K
-                      TEMP = ALPHA*B(J,L)
-                      DO 150 I = ISTART,ISTOP
-                          C(I,J) = C(I,J) + TEMP*A(I,L)
-  150                 CONTINUE
-  160             CONTINUE
-  170         CONTINUE
-          ELSE
-*
-*           Form  C := alpha*A**T*B**T + beta*C
-*
-              DO 200 J = 1,N
-                  IF (UPPER) THEN
-                      ISTART = 1
-                      ISTOP  = J
-                  ELSE
-                      ISTART = J
-                      ISTOP  = N
-                  END IF
-
-                  DO 190 I = ISTART, ISTOP
-                      TEMP = ZERO
-                      DO 180 L = 1,K
-                          TEMP = TEMP + A(L,I)*B(J,L)
-  180                 CONTINUE
-                      IF (BETA.EQ.ZERO) THEN
-                          C(I,J) = ALPHA*TEMP
-                      ELSE
-                          C(I,J) = ALPHA*TEMP + BETA*C(I,J)
-                      END IF
-  190             CONTINUE
-  200         CONTINUE
-          END IF
-      END IF
-*
-      RETURN
-*
-*     End of DGEMMTR
-*
-      END

BLAS/SRC/dgemv.f

@@ -155,7 +155,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dger.f

@@ -127,7 +127,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dnrm2.f90

@@ -85,12 +85,11 @@
 !> \endverbatim
 !>
 !  =====================================================================
-function DNRM2( n, x, incx )
-   implicit none
+function DNRM2( n, x, incx ) 
    integer, parameter :: wp = kind(1.d0)
    real(wp) :: DNRM2
 !
-!  -- Reference BLAS level1 routine --
+!  -- Reference BLAS level1 routine (version 3.9.1) --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021

BLAS/SRC/drot.f

@@ -89,7 +89,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/drotg.f90

@@ -89,7 +89,6 @@
 !
 !  =====================================================================
 subroutine DROTG( a, b, c, s )
-   implicit none
    integer, parameter :: wp = kind(1.d0)
 !
 !  -- Reference BLAS level1 routine --

BLAS/SRC/drotm.f

@@ -38,10 +38,6 @@
 *>    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:
@@ -97,7 +93,6 @@
 *
 *  =====================================================================
       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,    --

BLAS/SRC/drotmg.f

@@ -24,18 +24,14 @@
 *> \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.
+*>    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..
 *>
-*>    H HAS ONE OF THE FOLLOWING FORMS:
-*>
-*>      DFLAG=-1.D0     DFLAG=0.D0      DFLAG=1.D0      DFLAG=-2.D0
+*>    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.)
@@ -91,7 +87,6 @@
 *
 *  =====================================================================
       SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
@@ -200,7 +195,7 @@
                   DH11 = ONE
                   DH22 = ONE
                   DFLAG = -ONE
-               ELSE IF (DFLAG.EQ.ONE) THEN
+               ELSE
                   DH21 = -ONE
                   DH12 = ONE
                   DFLAG = -ONE
@@ -225,7 +220,7 @@
                   DH11 = ONE
                   DH22 = ONE
                   DFLAG = -ONE
-               ELSE IF (DFLAG.EQ.ONE) THEN
+               ELSE
                   DH21 = -ONE
                   DH12 = ONE
                   DFLAG = -ONE

BLAS/SRC/dsbmv.f

@@ -181,7 +181,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dscal.f

@@ -76,7 +76,6 @@
 *>
 *  =====================================================================
       SUBROUTINE DSCAL(N,DA,DX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/dsdot.f

@@ -116,7 +116,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dskewsymm.f

@@ -1,365 +0,0 @@
-*> \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

BLAS/SRC/dskewsymv.f

@@ -1,327 +0,0 @@
-*> \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

BLAS/SRC/dskewsyr2.f

@@ -1,294 +0,0 @@
-*> \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

BLAS/SRC/dskewsyr2k.f

@@ -1,395 +0,0 @@
-*> \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

BLAS/SRC/dspmv.f

@@ -144,7 +144,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dspr.f

@@ -124,7 +124,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dspr2.f

@@ -139,7 +139,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dswap.f

@@ -79,7 +79,6 @@
 *>
 *  =====================================================================
       SUBROUTINE DSWAP(N,DX,INCX,DY,INCY)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/dsymm.f

@@ -186,7 +186,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dsymv.f

@@ -149,7 +149,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dsyr.f

@@ -129,7 +129,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dsyr2.f

@@ -144,7 +144,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dsyr2k.f

@@ -189,7 +189,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dsyrk.f

@@ -166,7 +166,6 @@
 *>
 *  =====================================================================
       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,    --

BLAS/SRC/dtbmv.f

@@ -183,7 +183,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -198,6 +197,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -267,24 +270,28 @@
               KPLUS1 = K + 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
                       IF (J.GT.K) KX = KX + INCX
    40             CONTINUE
@@ -292,25 +299,29 @@
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
                       IF ((N-J).GE.K) KX = KX - INCX
    80             CONTINUE

BLAS/SRC/dtbsv.f

@@ -186,7 +186,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -201,6 +200,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -270,51 +273,59 @@
               KPLUS1 = K + 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
    20             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 40 J = N,1,-1
                       KX = KX - INCX
-                      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
+                      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
                       JX = JX - INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
                       KX = KX + INCX
-                      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
+                      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
                       JX = JX + INCX
    80             CONTINUE
               END IF

BLAS/SRC/dtpmv.f

@@ -139,7 +139,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -154,6 +153,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -216,25 +219,29 @@
               KK = 1
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
                       KK = KK + J
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
                       KK = KK + J
    40             CONTINUE
@@ -243,26 +250,30 @@
               KK = (N* (N+1))/2
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
                       KK = KK - (N-J+1)
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
                       KK = KK - (N-J+1)
    80             CONTINUE

BLAS/SRC/dtpsv.f

@@ -141,7 +141,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -156,6 +155,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -218,25 +221,29 @@
               KK = (N* (N+1))/2
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
                       KK = KK - J
    20             CONTINUE
               ELSE
                   JX = KX + (N-1)*INCX
                   DO 40 J = N,1,-1
-                      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
+                      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
                       JX = JX - INCX
                       KK = KK - J
    40             CONTINUE
@@ -245,25 +252,29 @@
               KK = 1
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
                       KK = KK + (N-J+1)
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
-                      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
+                      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
                       JX = JX + INCX
                       KK = KK + (N-J+1)
    80             CONTINUE

BLAS/SRC/dtrmm.f

@@ -174,7 +174,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -273,23 +272,27 @@
               IF (UPPER) THEN
                   DO 50 J = 1,N
                       DO 40 K = 1,M
-                          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
+                          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
    40                 CONTINUE
    50             CONTINUE
               ELSE
                   DO 80 J = 1,N
                       DO 70 K = M,1,-1
-                          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
+                          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
    70                 CONTINUE
    80             CONTINUE
               END IF
@@ -334,10 +337,12 @@
                           B(I,J) = TEMP*B(I,J)
   150                 CONTINUE
                       DO 170 K = 1,J - 1
-                          TEMP = ALPHA*A(K,J)
-                          DO 160 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  160                     CONTINUE
+                          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
   170                 CONTINUE
   180             CONTINUE
               ELSE
@@ -348,10 +353,12 @@
                           B(I,J) = TEMP*B(I,J)
   190                 CONTINUE
                       DO 210 K = J + 1,N
-                          TEMP = ALPHA*A(K,J)
-                          DO 200 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  200                     CONTINUE
+                          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
   210                 CONTINUE
   220             CONTINUE
               END IF
@@ -362,10 +369,12 @@
               IF (UPPER) THEN
                   DO 260 K = 1,N
                       DO 240 J = 1,K - 1
-                          TEMP = ALPHA*A(J,K)
-                          DO 230 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  230                     CONTINUE
+                          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
   240                 CONTINUE
                       TEMP = ALPHA
                       IF (NOUNIT) TEMP = TEMP*A(K,K)
@@ -378,10 +387,12 @@
               ELSE
                   DO 300 K = N,1,-1
                       DO 280 J = K + 1,N
-                          TEMP = ALPHA*A(J,K)
-                          DO 270 I = 1,M
-                              B(I,J) = B(I,J) + TEMP*B(I,K)
-  270                     CONTINUE
+                          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
   280                 CONTINUE
                       TEMP = ALPHA
                       IF (NOUNIT) TEMP = TEMP*A(K,K)

BLAS/SRC/dtrmv.f

@@ -144,7 +144,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -159,6 +158,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -225,45 +228,53 @@
           IF (LSAME(UPLO,'U')) THEN
               IF (INCX.EQ.1) THEN
                   DO 20 J = 1,N
-                      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)
+                      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
    20             CONTINUE
               ELSE
                   JX = KX
                   DO 40 J = 1,N
-                      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)
+                      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
                       JX = JX + INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = N,1,-1
-                      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)
+                      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
    60             CONTINUE
               ELSE
                   KX = KX + (N-1)*INCX
                   JX = KX
                   DO 80 J = N,1,-1
-                      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)
+                      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
                       JX = JX - INCX
    80             CONTINUE
               END IF

BLAS/SRC/dtrsm.f

@@ -178,7 +178,6 @@
 *>
 *  =====================================================================
       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,    --
@@ -211,8 +210,8 @@
       LOGICAL LSIDE,NOUNIT,UPPER
 *     ..
 *     .. Parameters ..
-      DOUBLE PRECISION ZERO
-      PARAMETER (ZERO=0.0D+0)
+      DOUBLE PRECISION ONE,ZERO
+      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 *     ..
 *
 *     Test the input parameters.
@@ -276,26 +275,34 @@
 *
               IF (UPPER) THEN
                   DO 60 J = 1,N
-                       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
+                      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
    50                 CONTINUE
    60             CONTINUE
               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
+                      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
    90                 CONTINUE
   100             CONTINUE
               END IF
@@ -334,33 +341,43 @@
 *
               IF (UPPER) THEN
                   DO 210 J = 1,N
-                      DO 170 I = 1,M
-                          B(I,J) = ALPHA*B(I,J)
-  170                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 170 I = 1,M
+                              B(I,J) = ALPHA*B(I,J)
+  170                     CONTINUE
+                      END IF
                       DO 190 K = 1,J - 1
-                          DO 180 I = 1,M
-                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  180                     CONTINUE
+                          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
   190                 CONTINUE
                       IF (NOUNIT) THEN
+                          TEMP = ONE/A(J,J)
                           DO 200 I = 1,M
-                              B(I,J) = B(I,J)/A(J,J)
+                              B(I,J) = TEMP*B(I,J)
   200                     CONTINUE
                       END IF
   210             CONTINUE
               ELSE
                   DO 260 J = N,1,-1
-                      DO 220 I = 1,M
-                          B(I,J) = ALPHA*B(I,J)
-  220                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 220 I = 1,M
+                              B(I,J) = ALPHA*B(I,J)
+  220                     CONTINUE
+                      END IF
                       DO 240 K = J + 1,N
-                          DO 230 I = 1,M
-                              B(I,J) = B(I,J) - A(K,J)*B(I,K)
-  230                     CONTINUE
+                          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
   240                 CONTINUE
                       IF (NOUNIT) THEN
+                          TEMP = ONE/A(J,J)
                           DO 250 I = 1,M
-                              B(I,J) = B(I,J)/A(J,J)
+                              B(I,J) = TEMP*B(I,J)
   250                     CONTINUE
                       END IF
   260             CONTINUE
@@ -372,34 +389,46 @@
               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) = B(I,K)/A(K,K)
+                              B(I,K) = TEMP*B(I,K)
   270                     CONTINUE
                       END IF
                       DO 290 J = 1,K - 1
-                          DO 280 I = 1,M
-                              B(I,J) = B(I,J) - A(J,K)*B(I,K)
-  280                     CONTINUE
+                          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
   290                 CONTINUE
-                      DO 300 I = 1,M
-                          B(I,K) = ALPHA*B(I,K)
-  300                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 300 I = 1,M
+                              B(I,K) = ALPHA*B(I,K)
+  300                     CONTINUE
+                      END IF
   310             CONTINUE
               ELSE
                   DO 360 K = 1,N
                       IF (NOUNIT) THEN
+                          TEMP = ONE/A(K,K)
                           DO 320 I = 1,M
-                              B(I,K) = B(I,K)/A(K,K)
+                              B(I,K) = TEMP*B(I,K)
   320                     CONTINUE
                       END IF
                       DO 340 J = K + 1,N
-                          DO 330 I = 1,M
-                              B(I,J) = B(I,J) - A(J,K)*B(I,K)
-  330                     CONTINUE
+                          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
   340                 CONTINUE
-                      DO 350 I = 1,M
-                          B(I,K) = ALPHA*B(I,K)
-  350                 CONTINUE
+                      IF (ALPHA.NE.ONE) THEN
+                          DO 350 I = 1,M
+                              B(I,K) = ALPHA*B(I,K)
+  350                     CONTINUE
+                      END IF
   360             CONTINUE
               END IF
           END IF

BLAS/SRC/dtrsv.f

@@ -140,7 +140,6 @@
 *
 *  =====================================================================
       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,    --
@@ -155,6 +154,10 @@
 *     ..
 *
 *  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION ZERO
+      PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
       DOUBLE PRECISION TEMP
@@ -221,44 +224,52 @@
           IF (LSAME(UPLO,'U')) THEN
               IF (INCX.EQ.1) THEN
                   DO 20 J = N,1,-1
-                      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
+                      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
    20             CONTINUE
               ELSE
                   JX = KX + (N-1)*INCX
                   DO 40 J = N,1,-1
-                      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
+                      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
                       JX = JX - INCX
    40             CONTINUE
               END IF
           ELSE
               IF (INCX.EQ.1) THEN
                   DO 60 J = 1,N
-                      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
+                      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
    60             CONTINUE
               ELSE
                   JX = KX
                   DO 80 J = 1,N
-                      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
+                      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
                       JX = JX + INCX
    80             CONTINUE
               END IF

BLAS/SRC/dzasum.f

@@ -69,7 +69,6 @@
 *>
 *  =====================================================================
       DOUBLE PRECISION FUNCTION DZASUM(N,ZX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/dznrm2.f90

@@ -86,12 +86,11 @@
 !> \endverbatim
 !>
 !  =====================================================================
-function DZNRM2( n, x, incx )
-   implicit none
+function DZNRM2( n, x, incx ) 
    integer, parameter :: wp = kind(1.d0)
    real(wp) :: DZNRM2
 !
-!  -- Reference BLAS level1 routine --
+!  -- Reference BLAS level1 routine (version 3.9.1) --
 !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 !     March 2021

BLAS/SRC/icamax.f90

@@ -1,193 +0,0 @@
-!> \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

BLAS/SRC/idamax.f

@@ -68,7 +68,6 @@
 *>
 *  =====================================================================
       INTEGER FUNCTION IDAMAX(N,DX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/isamax.f

@@ -68,7 +68,6 @@
 *>
 *  =====================================================================
       INTEGER FUNCTION ISAMAX(N,SX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/izamax.f90

@@ -1,193 +0,0 @@
-!> \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

BLAS/SRC/lsame.f

@@ -50,7 +50,6 @@
 *
 *  =====================================================================
       LOGICAL FUNCTION LSAME(CA,CB)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/sasum.f

@@ -69,7 +69,6 @@
 *>
 *  =====================================================================
       REAL FUNCTION SASUM(N,SX,INCX)
-      IMPLICIT NONE
 *
 *  -- Reference BLAS level1 routine --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --

BLAS/SRC/saxpby.f

@@ -1,148 +0,0 @@
-*> \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

BLAS/SRC/saxpy.f

@@ -86,7 +86,6 @@
 *>
 *  =====================================================================
       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,    --