diff --git a/.ci-pipelines/Dockerfile b/.ci-pipelines/Dockerfile deleted file mode 100644 index acbc2f7..0000000 --- a/.ci-pipelines/Dockerfile +++ /dev/null @@ -1,44 +0,0 @@ -# This dockerfile is used to run the CI unit tests on azure pipelines. -# To test build locally run the following commands: -# `docker build -f ./.ci-pipelines/Dockerfile --platform=amd64 -t kpp-build .` from top level of project -# To exec into the container run: `docker run -it --entrypoint /bin/bash kpp-build` - -# Set OS, a reference to maintainers, and the working directory. -# -# NOTE: Add --platform=amd64 (in .build-testing.yml and in your local -# aliases for 'docker build'. This will make sure that the C-I tests always -# run on AMD64 hardware. Otherwise, Docker will use the architecture of -# your current hardware, which might cause differences in library paths -# (e.g. /lib/aarch64-linux-gnu on Mac/M1 vs. /lib/x86_64-linux-gnu on Linux). -# -- Lucas Estrada & Bob Yantosca (06 Jul 2022) -FROM ubuntu:20.04 -LABEL maintainer="GEOS-Chem Support Team " -WORKDIR /kpp - -# Install KPP software dependencies -RUN apt-get update && apt-get install -y git gcc gfortran make flex bison - -COPY . . - -# Set environment variables for building KPP -ENV KPP_HOME=/kpp -ENV CC=gcc CXX=g++ FC=gfortran -ENV F90=$FC F9X=$FC F77=$FC -ENV KPP_FLEX_LIB_DIR=/lib/x86_64-linux-gnu - -# For C-I tests MAX_EQN and MAX_SPECIES must be < 1024: -# as this will fit within the memory limits on Azure. -RUN sed -i 's/#define MAX_EQN .*/#define MAX_EQN 1023/g' /kpp/src/gdata.h \ - && sed -i 's/#define MAX_SPECIES .*/#define MAX_SPECIES 1023/g' /kpp/src/gdata.h - -# Disable MCM test -RUN sed -i 's/DO_MCM=1/DO_MCM=0/g' /kpp/.ci-pipelines/ci-common-defs.sh - -# Build KPP executable and ensure testing scripts are executable -RUN cd /kpp/src/ && make -RUN chmod +x /kpp/.ci-pipelines/ci-common-defs.sh -RUN chmod +x /kpp/.ci-pipelines/ci-testing-script.sh - -# Run C-I tests -RUN /kpp/.ci-pipelines/ci-testing-script.sh - diff --git a/.ci-pipelines/build-testing.yml b/.ci-pipelines/build-testing.yml deleted file mode 100644 index 8d4c70d..0000000 --- a/.ci-pipelines/build-testing.yml +++ /dev/null @@ -1,27 +0,0 @@ -# Quick test pipeline: -# -# This pipeline runs several example mechanisms to check for errors -# impacting basic functionality of kpp - -# Run a C-I test when a push to any branch is made. -trigger: - branches: - include: - - '*' -pr: - branches: - include: - - '*' - -pool: - vmImage: 'ubuntu-24.04' -# Login to Docker Hub, build the image, and push the built image -# to Docker Hub -steps: - - task: Docker@2 - displayName: Build image - inputs: - command: build - buildContext: $(Build.Repository.LocalPath) # The path to the source code repo - Dockerfile: .ci-pipelines/Dockerfile - arguments: --platform=amd64 --progress=plain diff --git a/.ci-pipelines/ci-cleanup-script.sh b/.ci-pipelines/ci-cleanup-script.sh index 6e94669..fb9e346 100755 --- a/.ci-pipelines/ci-cleanup-script.sh +++ b/.ci-pipelines/ci-cleanup-script.sh @@ -17,8 +17,8 @@ cwd=$(pwd -P) for this_test in ${GENERAL_TESTS}; do clean_ci_test_folder "${this_test}" "${this_test}" "${cwd}" done -clean_ci_test_folder "${MCM_1}" "${MCM_1}" "${cwd}" -clean_ci_test_folder "${MCM_2}" "mcm" "${cwd}" +clean_ci_test_folder "${MCM_1}" mcm "${cwd}" +clean_ci_test_folder "${MCM_2}" mcm "${cwd}" # Remove any log files used to store C-I test results cd $cwd diff --git a/.ci-pipelines/ci-common-defs.sh b/.ci-pipelines/ci-common-defs.sh index 8f2c730..070861a 100755 --- a/.ci-pipelines/ci-common-defs.sh +++ b/.ci-pipelines/ci-common-defs.sh @@ -16,6 +16,7 @@ F90_feuler F90_lsode F90_radau F90_rk +F90_rkadj F90_rktlm F90_ros F90_rosadj @@ -27,7 +28,9 @@ F90_rostlm F90_ros_upcase F90_saprc_2006 F90_sd +F90_sd4 F90_sdadj +F90_sdtlm F90_seulex F90_small_strato " diff --git a/.github/workflows/lint-ci-workflows.yml b/.github/workflows/lint-ci-workflows.yml new file mode 100644 index 0000000..39a0f61 --- /dev/null +++ b/.github/workflows/lint-ci-workflows.yml @@ -0,0 +1,58 @@ +# Workflow to run linting checks on source +name: Lint + +# Controls when the workflow will run +on: + # Triggers the workflow on pushes to the "main" branch, i.e., PR merges + push: + branches: [ "main" ] + + # Triggers the workflow on pushes to open pull requests with code changes + pull_request: + paths: + - '.github/workflows/*.yml' + + # Allows you to run this workflow manually from the Actions tab + workflow_dispatch: + +# Cancel jobs running if new commits are pushed +concurrency: + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} + cancel-in-progress: true + +# Workflow run - one or more jobs that can run sequentially or in parallel +jobs: + # This workflow contains a single job called "lint" + lint: + # The type of runner that the job will run on + runs-on: ubuntu-latest + strategy: + fail-fast: false + + # Steps represent a sequence of tasks that will be executed as part of the job + steps: + # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it + - name: Checkout code + with: + persist-credentials: false + uses: actions/checkout@v6 + + - name: Install Python + uses: actions/setup-python@v6 + with: + python-version: '3.x' + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m venv ci_venv + . ci_venv/bin/activate + pip install zizmor==0.9.2 + + # Apply GitHub Actions linter, zizmor + - name: zizmor + if: always() + run: | + cd ${{ github.workspace }} + . ci_venv/bin/activate + zizmor .github/workflows/*.yml diff --git a/.github/workflows/run-ci-tests.yml b/.github/workflows/run-ci-tests.yml new file mode 100644 index 0000000..66fb966 --- /dev/null +++ b/.github/workflows/run-ci-tests.yml @@ -0,0 +1,68 @@ +--- +# +# GitHub Action to run C-I tests with GCC compilers +# +name: Run C-I tests + +on: + push: + branches: + - '**' + pull_request: + branches: + - '**' + +jobs: + build: + runs-on: ubuntu-latest + env: + KPP_HOME: ${{ github.workspace }} + KPP_FLEX_LIB_DIR: /usr/lib/x86_64-linux-gnu + strategy: + matrix: + gcc-version: [9, 10, 11, 12, 13] + + + name: Run C-I tests with GCC ${{ matrix.gcc-version }} + steps: + - name: Checkout KPP + uses: actions/checkout@v6 + with: + persist-credentials: false + + + - name: Install libraries + run: | + sudo apt update -y + sudo apt install -y flex bison libfl-dev + sudo apt install -y gcc-${{ matrix.gcc-version }} g++-${{ matrix.gcc-version }} gfortran-${{ matrix.gcc-version }} + sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-${{ matrix.gcc-version }} 100 \ + --slave /usr/bin/g++ g++ /usr/bin/g++-${{ matrix.gcc-version }} \ + --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${{ matrix.gcc-version }} + + + - name: Verify environment + run: | + echo '%%% COMPILERS %%%' + gcc --version + echo '' + g++ --version + echo '' + gfortran --version + echo '' + echo '%%% ENVIRONMENT %%%' + echo "flex = $(which flex)" + echo "bison = $(which bison)" + echo "KPP_HOME = $KPP_HOME" + echo "KPP_FLEX_LIB_DIR = $KPP_FLEX_LIB_DIR" + + - name: Compile with Make + run: | + cd $KPP_HOME/src + make all + + + - name: Run C-I tests + run: | + cd $KPP_HOME/.ci-pipelines + ./ci-testing-script.sh diff --git a/.gitignore b/.gitignore index c0ef5f1..dca92dc 100644 --- a/.gitignore +++ b/.gitignore @@ -1,73 +1,72 @@ +# KPP .gitignore file +# +# For info on how .gitignore rules work, see: # https://www.atlassian.com/git/tutorials/saving-changes/gitignore # https://linuxize.com/post/gitignore-ignoring-files-in-git -# temporary files: -*.~* -*~ -.\#* -\#* +# Keep empty directory for executables, but ignore its contents +/bin/* + +# Ignore temporary files everywhere +**/*.~* +**/*~ +**/.\#* +**/\#* -# files generated by bison: +# Ignore files generated by bison (in src/) /src/lex.yy.c /src/y.tab.c /src/y.tab.h -# from compiler: -*.mod -*.o -*.exe +# Ignore files generated by the compiler (everywhere) +**/*.mod +**/*.o +**/*.exe -# results from ci-tests: +# Ignore files generated by C-I tests (in ci-tests/*) +/ci-tests/**/*.dat +/ci-tests/**/*.h /ci-tests/**/*.m +/ci-tests/**/Makefile_* /ci-tests/**/*results* -/ci-tests/**/*.exe -/ci-tests/**/*.h -/ci-tests/**/*.log -/ci-tests/**/*.dat -# keep empty directory for executables, but not its contents: -/bin/* -!/bin/.gitkeep - -# KPP-generated files: -*_Function* -*_Global* -*_Hessian* -*_HessianSP* -*_Initialize* -*_Integrator* -*_Jacobian* -*_JacobianSP* -*_LinearAlgebra* -*.log -*_Main* -*_mex_Fun* -*_mex_Hessian* -*_mex_Jac_SP* -*_Model* -*_Monitor* -*_Parameters* -*_Precision* -*_Rates* -*_Stoichiom* -*_StoichiomSP* -*_Util* -Makefile_* +# Ignore files generated by KPP (everywhere) +**/*_Function* +**/*_Global* +**/*_Hessian* +**/*_HessianSP* +**/*_Initialize* +**/*_Integrator* +**/*_Jacobian* +**/*_JacobianSP* +**/*_LinearAlgebra* +**/*.log +**/*_Main* +**/*_*ex_Fun* +**/*_*ex_Hessian* +**/*_*ex_Jac_SP* +**/*_Model* +**/*_Monitor* +**/*_Parameters* +**/*_Precision* +**/*_Rates* +**/*_Stoichiom* +**/*_StoichiomSP* +**/*_Util* -# LaTeX and BibTeX files: -/manual/*.aux -/manual/*.bbl -/manual/*.blg -/manual/*.log -/manual/*.toc - -# ReadTheDocs files +# Ignore files generated by ReadTheDocs (in docs/build) docs/build/* -# are there any *.dat files to exclude? -*.dat - -# Other files/dirs to exclude +# Ignore other files & folders *.pdf -/examples/mcm*/__pycache__ +**/__pycache__ +# IMPORTANT: Do not ignore these files! +!/bin/.gitkeep +!/docs/figures/*.pdf +!/util/*.c +!/util/*.f +!/util/*.h +!/util/*.f90 +!/util/*.m +!/util/Makefile* diff --git a/CHANGELOG.md b/CHANGELOG.md index 62c6756..bf80d28 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,30 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [3.3.1] - 2026-03-27 +### Added +- Added GitHub Action to run C-I tests with GCC compilers v9, v10, v11, v12, and v13 +- Added "Lint" GitHub Action to check other actions for security issues +- Added new example files: `rkadj.kpp`, `sd4.kpp`, `sdtlm.kpp` +- Added new C-I tests: `F90_rkadj`, `F90_sd4`, `F90_sdtlm` +- Added `docs/read_the_docs_environment.yml` + +### Changed +- Updated ReadTheDocs documentation to reflect that C-I tests are now done as a GitHub Action +- Updated the `docs/source/reference/editing_these_docs.rst` instructions for clarity +- Updated GitHub Actions to the latest versions + +## Fixed +- Moved the `which kpp` instruction to the end of the "Build the KPP executable" section in the installation guide on ReadTheDocs +- Updated rules to ignore files in `.gitignore` and updated comments accordingly +- Fixed a bug that prevented `.ci-pipelines/ci-cleanup-script.sh` from removing KPP-generated files for MCM mechanisms +- Fixed typo in error message in `int/rosenbrock_autoreduce.f90` +- Fixed a bug in Makefile template that caused `KPP_ROOT_Stoichiom.f90` not compiled if `#STOICMAT ON` and `#HESSIAN OFF` + +### Removed +- Removed C-I tests on Microsoft Azure Dev Pipelines +- Replaced BLAS functions (`WAXPY`, `WCOPY`, `WSCAL`, `WADD`, `WLAMCH`, `WDOT`) with pure F90 code from `int/*.f90` integrators (thanks to AI for the help) + ## [3.3.0] - 2025-07-17 ### Added - New integrator: `rosenbrock_h211b_qssa.f90` diff --git a/ci-tests/F90_rkadj/F90_rkadj.kpp b/ci-tests/F90_rkadj/F90_rkadj.kpp new file mode 120000 index 0000000..b3bf812 --- /dev/null +++ b/ci-tests/F90_rkadj/F90_rkadj.kpp @@ -0,0 +1 @@ +../../examples/rkadj.kpp \ No newline at end of file diff --git a/ci-tests/F90_sd4/F90_sd4.kpp b/ci-tests/F90_sd4/F90_sd4.kpp new file mode 120000 index 0000000..f11975e --- /dev/null +++ b/ci-tests/F90_sd4/F90_sd4.kpp @@ -0,0 +1 @@ +../../examples/sd4.kpp \ No newline at end of file diff --git a/ci-tests/F90_sdtlm/F90_sdtlm.kpp b/ci-tests/F90_sdtlm/F90_sdtlm.kpp new file mode 120000 index 0000000..eb0b88c --- /dev/null +++ b/ci-tests/F90_sdtlm/F90_sdtlm.kpp @@ -0,0 +1 @@ +../../examples/sdtlm.kpp \ No newline at end of file diff --git a/docs/read_the_docs_environment.yml b/docs/read_the_docs_environment.yml new file mode 100644 index 0000000..5aefcac --- /dev/null +++ b/docs/read_the_docs_environment.yml @@ -0,0 +1,24 @@ +--- +# ====================================================================== +# ReadTheDocs environment file +# +# If you wish to build a Mamba/Conda environment with the dependencies +# needed for building the ReadTheDocs documentation, use: +# +# $ mamba env create -n rtd_env --file=read_the_docs_environment.yml +# ====================================================================== +name: rtd_env +channels: + - conda-forge + - nodefaults +dependencies: + - python==3.12.0 + - sphinx==7.2.6 + - sphinx_rtd_theme==2.0.0 + - sphinxcontrib-bibtex==2.6.2 + - sphinx-autobuild==2021.3.14 + - recommonmark==0.7.1 + - docutils==0.20.1 + - jinja2==3.1.6 + + diff --git a/docs/source/conf.py b/docs/source/conf.py index 4fbd6fd..411e75b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -25,7 +25,7 @@ # The full version, including alpha/beta/rc tags # (version numbers must be synchronized in CHANGELOG.md, src/gdata.h, # docs/source/conf.py and https://en.wikipedia.org/wiki/Kinetic_PreProcessor) -release = "3.3.0" +release = "3.3.1" # -- General configuration --------------------------------------------------- diff --git a/docs/source/getting_started/00_revision_history.rst b/docs/source/getting_started/00_revision_history.rst index 9591aaf..7868690 100644 --- a/docs/source/getting_started/00_revision_history.rst +++ b/docs/source/getting_started/00_revision_history.rst @@ -10,6 +10,22 @@ of the changes, read `CHANGELOG.md .. _unreleased: +.. _kpp331: + +========= +KPP 3.3.1 +========= + +- Replaced calls to BLAS functions in F90 integrators with core F90 + array operations +- Migrated C-I tests from Microsoft Azure Dev pipelines to GitHub + Actions +- Added new C-I tests: :literal:`F90_rkadj`, :literal:`F90_sd4`, + :literal:`F90_sdtlm` +- Fixed a compiliation issue that prevented + :literal:`KPP_ROOT_Stochiom.F90` from being compiled when + :literal:`#STOCHIOM ON` and :literal:`#HESSIAN OFF` are used + .. _kpp330: ========= diff --git a/docs/source/getting_started/01_installation.rst b/docs/source/getting_started/01_installation.rst index aa3ebf9..008180c 100644 --- a/docs/source/getting_started/01_installation.rst +++ b/docs/source/getting_started/01_installation.rst @@ -1,4 +1,4 @@ -.. _install: +.. _installation: ############ Installation @@ -6,227 +6,309 @@ Installation This section can be skipped if KPP is already installed on your system. +.. _installation-download: + ======================== -Download KPP from Github +Download KPP from GitHub ======================== -Clone the KPP source code from the `KPP Github repository +Clone the KPP source code from the `KPP GitHub repository `_: .. code-block:: console - $ cd $HOME - $ git clone https://github.com/KineticPreProcessor/KPP.git + cd $HOME + git clone https://github.com/KineticPreProcessor/KPP.git + +This will create a directory named :file:`KPP` in your home directory. -This will create a directory named KPP in your home directory. +.. _installation-env-vars: ======================================== Define the KPP_HOME environment variable ======================================== -Define the :envvar:`$KPP_HOME` environment variable to point to the -complete path where KPP is installed. Also, add the path of the KPP -executable to the :envvar:`$PATH` environment variable. +Define the :envvar:`KPP_HOME` environment variable to point to the +complete path where KPP was cloned. Also add the path of the KPP +executable to the :envvar:`PATH` environment variable. These commands +may be placed in your shell startup file, as shown below. -If you are using the Unix C-shell (aka :program:`csh`), add -add these statements to your :file:`$HOME/.cshrc` file: +.. important:: -.. code-block:: csh + The :file:`git clone` command above creates a directory named + :file:`KPP` (uppercase). Make sure that :envvar:`KPP_HOME` points + to this exact directory name. On case-sensitive Linux filesystems, + :file:`KPP` and :file:`kpp` are **different** directories. - setenv KPP_HOME $HOME/kpp - setenv PATH ${PATH}:$KPP_HOME/bin -and then apply the settings with: +**If you use bash:** -.. code-block:: console + Add these statements to your :file:`$HOME/.bashrc` file: - $ source $HOME/.cshrc + .. code-block:: bash -If, on the other hand, you are using the Unix :program:`bash` shell, -add these statements to your :file:`$HOME/.bashrc` file: + export KPP_HOME=$HOME/KPP + export PATH=$PATH:$KPP_HOME/bin -.. code-block:: bash + Then apply the settings with: - export KPP_HOME=$HOME/kpp - export PATH=$PATH:$KPP_HOME/bin + .. code-block:: console -and then apply the settings with: + $ source $HOME/.bashrc -.. code-block:: console +**If you use zsh (macOS Catalina and later)** - $ source $HOME/.bashrc + Add these statements to your :file:`$HOME/.zshrc` file: -Now if you type: + .. code-block:: zsh -.. code-block:: console + export KPP_HOME=$HOME/KPP + export PATH=$PATH:$KPP_HOME/bin + + Then apply the settings with: + + .. code-block:: console + + $ source $HOME/.zshrc + + + .. important:: + + macOS Catalina (2019) and later versions use :program:`zsh` + as the default shell instead of :program:`bash`. If you are + using :program:`bash` on macOS, source :file:`~/.bash_profile` + rather than :file:`~/.bashrc`, because macOS Terminal opens a + login shell by default and login shells do not source + :file:`~/.bashrc` automatically. The safest approach is to + add the following line to your :file:`~/.bash_profile`: - $ which kpp + .. code-block:: bash -the path to the executable file (:file:`kpp`) will be displayed. This -path should match the path specified by :file:`$KPP_HOME/bin/kpp`. + [ -f ~/.bashrc ] && source ~/.bashrc -.. _test-for-dependencies: +**If you use csh or tcsh:** -===================================================== -Test if KPP dependencies are installed on your system -===================================================== + Add these statements to your :file:`$HOME/.cshrc` file: -KPP depends on several other Unix packages. Before using KPP for the -first time, test if these are installed on your system. If any of -these packages are missing, you can install them with your -system's package manager (e.g. :program:`apt` for Ubuntu, -:program:`yum` for Fedora, :program:`homebrew` for MacOS, etc.), or -with `Spack `_. + .. code-block:: csh + + setenv KPP_HOME $HOME/KPP + setenv PATH ${PATH}:$KPP_HOME/bin + + Then apply the settings with: + + .. code-block:: console + + $ source $HOME/.cshrc + +.. _installation-dependencies: + +====================================== +Test if KPP dependencies are installed +====================================== + +KPP depends on several other Unix/Linux packages. Before using KPP +for the first time, test whether these are installed on your system. +If any packages are missing, you can install them with your system +package manager (e.g. :program:`apt` for Ubuntu, :program:`yum` for +Fedora, :program:`homebrew` for macOS) or with +`Spack `_. + +.. _installation-gcc: gcc --- .. important:: - You might have to follow some :ref:`additional configuration - and installation steps ` regarding - :program:`gcc` on MacOS X systems. + macOS users: please read :ref:`installation-macos-gcc` before + proceeding. -KPP uses the `GNU Compiler Collection `_ (aka -:program:`gcc`) by default. A version of :program:`gcc` comes -pre-installed with most Linux or MacOS systems. To test if -:program:`gcc` is installed on your system, type: +KPP uses the `GNU Compiler Collection `_ +(:program:`gcc`) by default. A version of :program:`gcc` comes +pre-installed with most Linux systems. To test whether +:program:`gcc` is installed, type: -.. code-block :: console +.. code-block:: console - $ gcc --version + gcc --version -This will display the version information, such as: +This will display version information such as: -.. code-block:: console +.. code-block:: none gcc (GCC) 11.2.0 Copyright (C) 2021 Free Software Foundation, Inc. - This is free software; see the source for copying conditions. There is NO - warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + This is free software; see the source for copying conditions. + There is NO warranty; not even for MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. + +.. _installation-sed: sed --- -The :program:`sed` utility is used to search for and replace text -in files. To test if :program:`sed` has been installed, type: + +The :program:`sed` utility searches for and replaces text in files. +To test whether :program:`sed` is installed, type: .. code-block:: console - $ which sed + which sed This will print the path to :program:`sed` on your system. +.. _installation-bison: + bison ----- -The :program:`bison` utility parses the chemical mechanism file into a -computer-readable syntax. To test :program:`bison` is installed, type: +The :program:`bison` utility parses the chemical mechanism file into +a computer-readable syntax. To test whether :program:`bison` is +installed, type: .. code-block:: console - $ which bison + which bison This will print the path to :program:`bison` on your system. -.. _flex-dep: +.. _installation-flex: flex ---- .. important:: - You might have to follow some :ref:`additional configuration - and installation steps ` regarding - :program:`flex` on MacOS X systems. + macOS users: please read :ref:`installation-macos-flex` before + proceeding. -The :program:`flex` (the Fast Lexical Analyzer) creates a scanner that -can recognize the syntax generated by :program:`bison`. To test if -:program:`flex` is installed, type: +The :program:`flex` (Fast Lexical Analyzer) utility creates a scanner +that recognizes the syntax generated by :program:`bison`. To test +whether :program:`flex` is installed, type: .. code-block:: console - $ which flex + which flex This will print the path to :program:`flex` on your system. You will also need to specify the path to the :program:`flex` library -files (:file:`libfl.so` or :file:`libfl.a`) in order to :ref:`build -the KPP executable `. This can be done with the -:program:`find` command: +files (:file:`libfl.so` or :file:`libfl.a`) in order to build the KPP +executable. Use the :program:`find` command to locate them: .. code-block:: console - $ find /usr/ -name "*libfl*" -print + find /usr/ -name "*libfl*" -print -This will generate a list of file paths such as shown below. Look for -the text :file:`libfl.`: +This generates output similar to the following. Look for the entry +containing :file:`libfl.`: -.. code-block:: console +.. code-block:: none /usr/include/libflashrom.h /usr/lib/gthumb/extensions/libflicker.so /usr/lib/gthumb/extensions/libflicker_utils.so /usr/lib/libflashrom.so.1.0.0 /usr/lib/libfl.so # <---- This is the flex library file - # ... etc ... + ... etc ... -Once you have located the directory where flex library file -resides (which in this example is :file:`/usr/lib`), use it to define -the :envvar:`KPP_FLEX_LIB_DIR` environment variable in your -:file:`.bashrc` (or :file:`.bash_aliases` file if you have one): +Once you have located the directory containing the flex library file +(in this example :file:`/usr/lib`), define the +:envvar:`KPP_FLEX_LIB_DIR` environment variable in your shell startup +file as shown below. -.. code-block:: bash +**If you use bash:** - export KPP_FLEX_LIB_DIR=/usr/lib - export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}:" + Add this code to your :file:`~/.bashrc` file: -Then apply the changes with: + .. code-block:: bash -.. code-block:: console + export KPP_FLEX_LIB_DIR=/usr/lib + export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}" + + Then apply the changes with: + + .. code-block:: console + + $ source ~/.bashrc + +**If you use zsh:** - . ~/.bashrc + Add this code to your :file:`~/.zshrc` file: + + .. code-block:: bash + + export KPP_FLEX_LIB_DIR=/usr/lib + export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}" + + Then apply the changes with: + + .. code-block:: console -KPP will use the path specified by :envvar:`KPP_FLEX_LIB_DIR` during -the compilation sequence (described in the next section). + $ source ~/.zshrc -.. _build-kpp-exec: +**If you use csh or tcsh:** + + Add this code to your :file:`~/.cshrc` file: + + .. code-block:: csh + + setenv KPP_FLEX_LIB_DIR /usr/lib + setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR} + + Then apply the changes with: + + .. code-block:: console + + $ source ~/.cshrc + +KPP will use the path specified by :envvar:`KPP_FLEX_LIB_DIR` +in the following section. + +.. _installation-build: ======================== Build the KPP executable ======================== -Change to the KPP/src directory: +Change to the :file:`KPP/src` directory: .. code-block:: console - $ cd $KPP_HOME/src + cd $KPP_HOME/src + +.. note:: -To clean a previously-built KPP installation, delete the KPP object -files and all the examples with: + The following :program:`make clean` and :program:`make distclean` + commands are only necessary if you have previously built KPP and + wish to start from a clean state. You can skip them on a fresh + clone. + +To remove previously-built KPP object files and example output: .. code-block:: console - $ make clean + make clean -To delete a previoulsy-built KPP executable as well, type: +To also remove a previously-built KPP executable: .. code-block:: console - $ make distclean + make distclean -KPP will use :program:`gcc` as the default compiler. If you would -like to use a different compiler (e.g. :program:`icc`), then edit -:file:`src/Makefile.defs` to add your compiler name. +KPP uses :program:`gcc` as the default compiler. If you would like +to use a different compiler (e.g. :program:`icc`), edit +:file:`src/Makefile.defs` to specify your compiler. Create the KPP executable with: .. code-block:: console - $ make + make You should see output similar to: -.. code-block:: console +.. code-block:: none gcc -g -Wall -Wno-unused-function -I/usr/include -c code.c gcc -g -Wall -Wno-unused-function -I/usr/include -c code_c.c @@ -242,214 +324,321 @@ You should see output similar to: gcc -g -Wall -Wno-unused-function -I/usr/include -c scanner.c gcc -g -Wall -Wno-unused-function -I/usr/include -c scanutil.c gcc -g -Wall -Wno-unused-function -I/usr/include -c y.tab.c - gcc -g -Wall -Wno-unused-function code.o code_c.o - code_f77.o code_f90.o code_matlab.o debug.o gen.o kpp.o + gcc -g -Wall -Wno-unused-function code.o code_c.o \ + code_f77.o code_f90.o code_matlab.o debug.o gen.o kpp.o \ lex.yy.o scanner.o scanutil.o y.tab.o -L/usr/lib -lfl -o kpp -This will create the executable file :file:`$KPP_HOME/bin/kpp`. +This creates the executable file :file:`$KPP_HOME/bin/kpp`. + +To verify that the KPP executable is accessible, type: -.. _additional-steps-macosx: +.. code-block:: console -============================== -Instructions for MacOS X users -============================== + which kpp -When installing KPP on a MacOS X system, some additional configuration -and installation steps may be necessary. +The path displayed should match :file:`$KPP_HOME/bin/kpp`. -.. _force-macos-to-recognize-gcc-compiler: +.. note:: -Force MacOS to recognize the gcc compiler ------------------------------------------ + This check must be performed **after** running :command:`make`. + The executable does not exist before the build step completes. -On MacOS X, if you type: +.. _installation-macos: -.. code-block:: console +========================== +Additional steps for macOS +========================== + +When installing KPP on a macOS system, some additional configuration +steps are necessary. The sections below address each one. + +.. _installation-macos-gcc: - $ gcc --version +Force macOS to use the GNU gcc compiler +--------------------------------------- -you will probably see output similar to: +On macOS, the command: .. code-block:: console + gcc --version + +will typically produce output similar to: + +.. code-block:: none + Apple clang version 13.1.6 (clang-1316.0.21.2.5) Target: x86_64-apple-darwin21.5.0 Thread model: posix InstalledDir: /Library/Developer/CommandLineTools/usr/bin -This is because MacOS X installs :program:`clang` as :program:`gcc`. -To force MacOS X to recognize the :program:`gcc` compiler, follow -these steps: +This is because macOS installs :program:`clang` as :program:`gcc`. +Follow the steps below to make your shell use the GNU +:program:`gcc` compiler instead. -#. Use the :program:`homebrew` package manager to install - :program:`gcc`: +1. Install gcc with Homebrew +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - .. code-block:: console +.. code-block:: console - $ brew install gcc + brew install gcc -#. Type this command: +2. Determine which version of gcc was installed +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - .. code-block:: console +.. code-block:: console - $ ls /usr/local/Cellar/gcc/*/bin/ | grep gcc + ls $(brew --prefix gcc)/bin/ | grep "^gcc-" - You should see output such as: +You will see output such as: - .. code-block:: console +.. code-block:: none + + gcc-13 + gcc-ar-13 + gcc-nm-13 + gcc-ranlib-13 + +This indicates that gcc version 13 was installed and that the +executable is named :file:`gcc-13`. Your installed version may differ. + +.. important:: + + The version number returned by the command above may differ from + the examples shown here. Always use the actual version number + found on your system in the steps that follow. - gcc-11* - gcc-ar-11* - gcc-nm-11* - gcc-ranlib-11* - # ... etc ... +3. Define compiler settings in your shell startup file +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - This output indicates :program:`gcc` major version 11 has been - installed, and that the gcc executable is called :code:`gcc-11`. - (Your version may differ.) +The following code block uses :command:`brew --prefix` to locate the +Homebrew installation directory automatically. This works correctly +on both Intel Macs (:file:`/usr/local`) and Apple Silicon Macs +(:file:`/opt/homebrew`). -#. Add the following code block to your :file:`.bashrc` file (or to your - :file:`.bash_aliases` file if you have one). This will define - aliases that will override :program:`clang` with :program:`gcc`. +**If you use bash:** + + Add the following to your :file:`~/.bashrc` (or + :file:`~/.bash_profile` if you do not have a + :file:`~/.bashrc`): .. code-block:: bash - #============================================================================ - # Compiler settings (MacOS) + #================================================================== + # Compiler settings (macOS) # - # NOTE: MacOSX installs Clang as /usr/bin/gcc, so we have to manually - # force reference to gcc-11, g++-11, and gfortran-11, which HomeBrew - # installs to /usr/local/bin. (bmy, 10/28/21) - #============================================================================ - alias gcc=gcc-11 - alias g++=g++-11 - alias gfortran=gfortran-11 - export CC=gcc - export CXX=g++-11 - export FC=gfortran-11 - export F77=gfortran-11 + # Detect the gcc version installed by Homebrew and define + # aliases and environment variables accordingly. + # This overrides the Apple clang that macOS installs as gcc. + #================================================================== + _GCC_VER=$(ls "$(brew --prefix gcc)/bin/" \ + | grep '^gcc-[0-9]' \ + | grep -oP '[0-9]+$' \ + | sort -n | tail -1) + + if [[ -n "${_GCC_VER}" ]]; then + alias gcc="gcc-${_GCC_VER}" + alias g++="g++-${_GCC_VER}" + alias gfortran="gfortran-${_GCC_VER}" + export CC="gcc-${_GCC_VER}" + export CXX="g++-${_GCC_VER}" + export FC="gfortran-${_GCC_VER}" + export F77="gfortran-${_GCC_VER}" + fi + unset _GCC_VER Then apply the changes with: .. code-block:: console - $ . ~/.bashrc + $ source ~/.bashrc -#. To check if your shell now recognizes the :program:`gcc` compiler, type: +**If you use zsh:** - .. code-block:: console + Add the following to your :file:`~/.zshrc`: - $ gcc --version + .. code-block:: zsh - You should see output similar to: + #================================================================== + # Compiler settings (macOS) + # + # Detect the gcc version installed by Homebrew and define + # aliases and environment variables accordingly. + # This overrides the Apple clang that macOS installs as gcc. + #================================================================== + _GCC_VER=$(ls "$(brew --prefix gcc)/bin/" \ + | grep '^gcc-[0-9]' \ + | grep -oE '[0-9]+$' \ + | sort -n | tail -1) + + if [[ -n "${_GCC_VER}" ]]; then + alias gcc="gcc-${_GCC_VER}" + alias g++="g++-${_GCC_VER}" + alias gfortran="gfortran-${_GCC_VER}" + export CC="gcc-${_GCC_VER}" + export CXX="g++-${_GCC_VER}" + export FC="gfortran-${_GCC_VER}" + export F77="gfortran-${_GCC_VER}" + fi + unset _GCC_VER + + Then apply the changes with: .. code-block:: console - gcc-11 (Homebrew GCC 11.3.0_1) 11.3.0 - Copyright (C) 2021 Free Software Foundation, Inc. - This is free software; see the source for copying conditions. There is NO - warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + $ source ~/.zshrc - This now indicates that your compiler is :program:`gcc` and not - :program:`clang`. +4. Verify that your shell now recognizes the GNU gcc compiler +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. _install-flex-with-homebrew: +.. code-block:: console -Install flex with homebrew --------------------------- + gcc --version + +You should see output similar to: -If your MacOS X computer does not have the :program:`flex` library -installed, then you can install it with :program:`homebrew`: +.. code-block:: none -.. code-block:: console + gcc-13 (Homebrew GCC 13.2.0) 13.2.0 + Copyright (C) 2023 Free Software Foundation, Inc. + This is free software; see the source for copying conditions. + There is NO warranty; not even for MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. - $ brew install flex +This confirms that your shell is using GNU :program:`gcc` rather +than :program:`clang`. -Unlike Linux pacakge managers, which would install the :program:`flex` -library files in the path :file:`/usr/lib/`, -:program:`homebrew` will install it to a path such as -:file:`/usr/local/Cellar/flex/X.Y.Z/lib/`. +.. _installation-macos-flex: -To find the version of :program:`flex` that has been installed by -:program:`homebrew`, type: +Install and configure flex on macOS +------------------------------------ + +If your macOS system does not have the :program:`flex` library +installed, install it with Homebrew: .. code-block:: console - $ ls /usr/local/Cellar/flex + brew install flex -and you will get a listing such as: +Unlike Linux package managers, which install the flex library files +under :file:`/usr/lib/`, Homebrew installs them under its own prefix +directory. -.. code-block:: console +Define the :envvar:`KPP_FLEX_LIB_DIR` environment variable using +:command:`brew --prefix` so that the path is determined automatically +regardless of the Homebrew prefix or the installed flex version: - 2.6.4_2 +**If you use bash:** -This indicates that the version of :program:`flex` on your system is -:code:`2.6.4_2` (the :code:`_2` denotes the number of bug-fix updates -since version :code:`2.6.4` was released). + Add the following to your :file:`~/.bashrc` (or + :file:`~/.bash_profile`): -The :program:`flex` library files (:file:`libfl.so` or -:file:`libfl.a`) will be found in :file:`lib/` subfolder. In this -example, the path will be: + .. code-block:: bash -.. code-block:: console + export KPP_FLEX_LIB_DIR=$(brew --prefix flex)/lib + export DYLD_LIBRARY_PATH="${DYLD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}" + + Then apply the changes with: + + .. code-block:: console + + $ source ~/.bashrc - /usr/local/Cellar/flex/2.6.4_2/lib +**If you use zsh:** -Knowing this, you can now define the :envvar:`KPP_FLEX_LIB_DIR` -environment variable :ref:`as described above `: + Add the following to your :file:`~/.zshrc`: -.. code-block:: bash + .. code-block:: zsh - export FLEX_LIB_DIR=/usr/local/Cellar/flex/2.6.4_2/lib + export KPP_FLEX_LIB_DIR=$(brew --prefix flex)/lib + export DYLD_LIBRARY_PATH="${DYLD_LIBRARY_PATH}:${KPP_FLEX_LIB_DIR}" + + Then apply the changes with: + + .. code-block:: console + + $ source ~/.zshrc + +.. note:: + + macOS uses :envvar:`DYLD_LIBRARY_PATH` to locate shared libraries + at runtime, not :envvar:`LD_LIBRARY_PATH` (which is the Linux + equivalent). If you find that :envvar:`DYLD_LIBRARY_PATH` is + blocked by System Integrity Protection (SIP), you can instead force + a hard link: + + .. code-block:: console -.. _macosx-limited-stack: + brew link --force flex + +.. _installation-macos-stack: Request maximum stack memory ----------------------------- +----------------------------- + +macOS imposes a hard limit on stack memory of 65532 bytes, which is +considerably less than what is available on typical GNU/Linux systems. +To ensure that KPP uses the maximum available stack memory, add the +following line to your shell startup file: -MacOS X has a hard limit of 65332 bytes for stack memory. This is -much less memory than what is available on GNU/Linux operating systems -such as Ubuntu, Fedora, etc. +**If you use bash:** -To make sure you are using the maximum amount of stack memory on MacOS -X add this command to your :file:`.bashrc` file: + Add this code to your :file:`~/.bashrc` file: -.. code-block:: bash + .. code-block:: bash - ulimit -s 65532 + ulimit -s 65532 -and then apply the change with: + Then apply the change with: -.. code-block:: console + .. code-block:: console + + $ source ~/.bashrc + +**If you use zsh:** + + Add this code to your :file:`~/.zshrc` file: - $ . ~/.bashrc + .. code-block:: zsh -This stack memory limit means that KPP will not be able to parse -mechanisms with more than about 2000 equations and 1000 species. -Because of this, we have added an :code:`#ifdef` block to KPP header -file :file:`src/gdata.h` to define the :code:`MAX_EQN` and -:code:`MAX_SPECIES` parameters accordingly: + ulimit -s 65532 -.. code-block:: C + Then apply the change with: + + .. code-block:: console + + $ source ~/.zshrc + +This stack memory restriction means that KPP cannot parse mechanisms +with more than approximately 2000 equations and 1000 species. To +account for this, the KPP header file :file:`src/gdata.h` defines the +:c:macro:`MAX_EQN` and :c:macro:`MAX_SPECIES` parameters conditionally: + +.. code-block:: c #ifdef MACOS - #define MAX_EQN 2000 // Max number of equations (MacOS only) - #define MAX_SPECIES 1000 // Max number of species (MacOS only) + #define MAX_EQN 2000 // Max number of equations (macOS only) + #define MAX_SPECIES 1000 // Max number of species (macOS only) #else - #define MAX_EQN 11000 // Max number of equations - #define MAX_SPECIES 6000 // Max number of species + #define MAX_EQN 11000 // Max number of equations + #define MAX_SPECIES 6000 // Max number of species #endif -If you find that KPP will not parse your mechanism, you can increase -:code:`MAX_EQN` and decrease :code:`MAX_SPECIES` (or vice-versa) as -needed, and then :ref:`rebuild the KPP executable `. +If KPP cannot parse your mechanism you can adjust :c:macro:`MAX_EQN` +and :c:macro:`MAX_SPECIES` in :file:`src/gdata.h` and then rebuild +the KPP executable. + +.. _installation-macos-case: -.. _macosx-case-insensitive: +macOS is case-insensitive +-------------------------- -Know that MacOS X is case-insenstive -------------------------------------- +.. warning:: -If you have two files with identical names except for case -(e.g. :file:`integrator.F90` and :file:`integrator.f90`) then MacOS X -will not be able to tell them apart. Because of this, you may -encounter an error if you try to commit such files into Git, etc. + macOS uses a case-insensitive filesystem by default. If your + project contains two files whose names differ only in case (for + example :file:`integrator.F90` and :file:`integrator.f90`), macOS + will treat them as the same file. This can cause build failures + and unexpected behaviour when working with version control systems + such as Git. Ensure that no two files in your project share a + case-insensitive name. diff --git a/docs/source/getting_started/02_running_kpp_sample_mech.rst b/docs/source/getting_started/02_running_kpp_sample_mech.rst index c3a0cb0..afdb8fc 100644 --- a/docs/source/getting_started/02_running_kpp_sample_mech.rst +++ b/docs/source/getting_started/02_running_kpp_sample_mech.rst @@ -261,7 +261,7 @@ You should see output similar to: .. code-block:: console - This is KPP-3.0.0. + This is KPP-X.Y.Z. KPP is parsing the equation file. KPP is computing Jacobian sparsity structure. diff --git a/docs/source/index.rst b/docs/source/index.rst index b91e1b6..1133c2b 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -31,9 +31,8 @@ The Kinetic PreProcessor: KPP This site provides instructions for :program:`KPP`, the Kinetic PreProcessor. Contributions (e.g., suggestions, edits, revisions) would be greatly -appreciated. See :ref:`editing-this-user-guide` and our -contributing guidelines. If you find something hard to -understand---let us know! +appreciated. See :ref:`editing_this_user_guide` and our contributing +guidelines. If you find something hard to understand---let us know! .. toctree:: :caption: Getting Started diff --git a/docs/source/reference/editing_these_docs.rst b/docs/source/reference/editing_these_docs.rst index 38df523..1f8df05 100644 --- a/docs/source/reference/editing_these_docs.rst +++ b/docs/source/reference/editing_these_docs.rst @@ -1,132 +1,372 @@ +.. |br| raw:: html -.. _editing-this-user-guide: +
+ +.. _editing_this_user_guide: ####################### Editing this User Guide ####################### -This user guide is generated with `Sphinx `_. -Sphinx is an open-source Python project designed to make writing software documentation easier. -The documentation is written in a reStructuredText (it's similar to markdown), which Sphinx extends for software documentation. -The source for the documentation is the :file:`docs/source` directory in top-level of the source code. +This user guide is generated with `Sphinx +`_. Sphinx is an open-source Python +project designed to make writing software documentation easier. The +documentation is written in :ref:`reStructuredText (reST) +`, a plaintext markup language that +Sphinx extends for software documentation. The source for the +documentation is the :file:`docs/source` directory in top-level of the +source code (and its subdirectories). + +.. _editing_this_user_guide_quickstart: =========== Quick start =========== -To build this user guide on your local machine, you need to install Sphinx. Sphinx is a Python 3 package and -it is available via :program:`pip`. This user guide uses the Read The Docs theme, so you will also need to -install :literal:`sphinx-rtd-theme`. It also uses the `sphinxcontrib-bibtex `_ -and `recommonmark `_ extensions, which you'll need to install. +First-time setup: Install Sphinx +-------------------------------- + +To build this user guide on your local machine, you need to install +Sphinx and its dependencies, which are listed in the table below. + +.. list-table:: + :header-rows: 1 + :widths: 30 50 20 + + * - Package + - Description + - Version + * - sphinx + - Creates online user manual documentation from markup text files + - 7.2.6 + * - `sphinx-autobuild `_ + - Dynamically builds Sphinx documentation and displays it in a + browser + - 2021.3.14 + * - `sphinx_rtd_theme `_ + - Sphinx theme for ReadTheDocs + - 2.0.0 + * - `sphinxcontrib-bibtex `_ + - Inserts LaTeX-style bibliography citations into ReadTheDocs + documentation + - 2.6.2 + * - `docutils `_ + - Processes plaintext documentation into HTML and other formats + - 0.20.1 + * - `recommonmark `_ + - Parses text for docutils + - 0.7.1 + * - `jinja2 `_ + - Replaces tokenized strings with text + - 3.1.6 + +We recommend that you create a standalone :program:`Conda` environment +to install Sphinx and its dependencies. The YAML file +:file:`docs/read_the_docs_environment.yml` contains the proper package +specifications. Use these commands: .. code-block:: console - $ cd $KPP_HOME/docs - $ pip install -r requirements.txt + $ cd docs + $ conda env create -n rtd_env --file=read_the_docs_environment.yml -Installing with :file:`requirements.txt` will make sure that the -proper versions of Sphinx and its dependencies will be installed. +This step only needs to be done once. -To build this user guide locally, navigate to the :file:`docs/` directory and make the :literal:`html` target. +Build the documentation +----------------------- -.. code-block:: console +#. Activate the :program:`Conda` environment containing + :program:`Sphinx` and its dependencies: + + .. code-block:: console + + $ conda activate rtd_env + +#. Navigate to the :file:`docs/` folder: + + .. code-block:: console + + (rtd_env) $ cd docs # Skip if you are already in the docs folder + +#. Check out the :file:`docs/dev` branch of this repository, as this + is the branch from which the :program:`latest` ReadTheDocs version + will be built: + + .. code-block:: console + + (rtd_env) $ git checkout docs/dev # Skip if you are already on the docs/dev branch + +#. Start the :command:`sphinx-autobuild` server: + + .. code-block:: console + + (rtd_env) $ sphinx-autobuild source build/html + + This will parse the reST-format files in the :file:`docs/source/` + directory tree and generate new HTML files in + :file:`docs/build/html`. |br| + |br| + +#. Remove any HTML files (in :file:`docs/build/html`) that might be + left behind from a previous build: - $ cd $KPP_HOME/docs - $ make html + .. code-block:: console + (rtd_env) $ make clean -This will build the user guide in :file:`docs/build/html`, and you can -open :file:`index.html` in your web-browser. The source files for the -user guide are found in :file:`docs/source`. +#. Open a web browser and navigate to :file:`localhost:8000`. |br| + |br| -.. note:: +#. Open your favorite text editor and start making changes to the + reST-format documentation files in the :file:`docs/source` + directory tree. While :program:`sphinx-autobuild` is running, you + will see your updates rendered in the web browser as soon as you + soon as you save your changes to disk. |br| + |br| - You can clean the documentation with :code:`make clean`. +#. Once you are satisfied with your edits, commit your changes to Git + and push the documentation to the :file:`docs/dev` remote branch of + this repository, |br| + |br| -============= -Learning reST -============= +#. Remove the generated HTML documentation files: -Writing reST can be tricky at first. Whitespace matters, and some directives -can be easily miswritten. Two important things you should know right away are: + .. code-block:: console -* Indents are 3-spaces -* "Things" are separated by 1 blank line. For example, a list or + (rtd_env) $ make clean + +#. Halt the :program:`sphinx-autobuild` server by typing + :program:`CTRL-C`. |br| + |br| + +#. Deactivate the :program:`Conda` environment: + + .. code-block:: console + + (rtd_env) $ conda deactivate + +.. _editing_this_user_guide_rest: + +========================== +Learning reStructured Text +========================== + +ReadTheDocs documentation is generated from text files in **reStructured +Text (reST)**, which is an easy-to-read, what-you-see-is-what-you-get +plaintext markup language. It is the default markup language used by +Sphinx. + +Writing reST can be tricky at first. Whitespace matters, and some +directives can be easily miswritten. Two important things you should +know right away are: + +- Indents are 3-spaces +- "Things" are separated by 1 blank line. For example, a list or code-block following a paragraph should be separated from the paragraph by 1 blank line. You should keep these in mind when you're first getting started. Dedicating an hour to learning reST will save you time in the -long-run. Below are some good resources for learning reST. +long-run. Below are some good resources for learning reST. -* `reStructuredText primer +- `reStructuredText primer `_: - (single best resource; however, it's better read than skimmed) -* Official `reStructuredText reference + (single best resource; however, it's better read than skimmed) |br| + |br| + +- Official `reStructuredText reference `_ - (there is *a lot* of information here) -* `Presentation by Eric Holscher + (there is *a lot* of information here) |br| + |br| + +- `Presentation by Eric Holscher `_ (co-founder of Read The Docs) at DjangoCon US 2015 (the entire presentation is good, but - reST is described from 9:03 to 21:04) -* `YouTube tutorial by Audrey Tavares's + reST is described from 9:03 to 21:04) |br| + |br| + +- `YouTube tutorial by Audrey Tavares `_ A good starting point would be Eric Holscher's presentations followed by the reStructuredText primer. +.. _editing_this_user_guide_style: + ================ Style guidelines ================ -.. important:: +This user guide is written in semantic markup. This is important so +that the user guide remains maintainable. Before contributing to +this documentation, please review our style guidelines +(below). When editing the source, please refrain from using +elements with the wrong semantic meaning for aesthetic +reasons. Aesthetic issues can be addressed by changes to the theme. + +Titles and headers +------------------ + +.. list-table:: + :header-rows: 1 + :widths: 40 60 + + * - Element + - reST Markup + * - Section header |br| (aka "Heading 1) + - Overline by :literal:`#` and underline by :literal:`#` + * - Sub-section header |br| (aka "Heading 2") + - Overline by :literal:`=` and underline by :literal:`=` + * - Sub-sub-section header |br| (aka "Heading 3") + - Underline by :literal:`-` + * - Sub-sub-sub-section header |br| (aka "Heading 4") + - Underline by :literal:`~` + * - Sub-sub-sub-sub-section header |br| (aka "Heading 5") + - Underline by :literal:`^` + +References and links +-------------------- + +.. list-table:: + :header-rows: 1 + :widths: 30 35 35 + + * - Element + - reST Markup Example + - Rendered text + * - Reference to a named anchor + - ``:ref:`editing_this_user_guide_quickstart``` + - :ref:`editing_this_user_guide_quickstart` + * - Renamed reference to a named anchor + - ``:ref:`Getting Started `` + - :ref:`Getting Started ` + * - HTML link + - ```ReadTheDocs `_`` + - `GEOS-Chem Manual `_ + +Other common style elements +--------------------------- + +.. list-table:: + :header-rows: 1 + :widths: 30 35 35 + + * - Element + - reST Markup Example + - Rendered text + * - File paths + - ``:file:`myfile.nc``` + - :file:`myfile.nc` + * - Directories + - ``:file:`/usr/bin``` + - :file:`/usr/bin` + * - Program names + - ``:program:`cmake``` + - :program:`cmake` + * - OS-level commands + - ``:program:`rm -rf``` + - :program:`rm -rf` + * - Environment variables + - ``:envvar:`$HOME``` + - :envvar:`$HOME` + * - Inline code or code variables + - ``:code:`PRINT*, "HELLO!"``` + - :code:`PRINT*, "HELLO!"` + * - Inline literal text + - ``:literal:`$``` + - :literal:`$` + +Indented code and text blocks +----------------------------- + +Code snippets should use :literal:`.. code-block:: ` +directives: + +Python +~~~~~~ + +.. code-block:: none + + .. code-block:: python + + import gcpy + print("hello world") + +Renders as: - This user guide is written in semantic markup. This is important so - that the user guide remains maintainable. Before contributing to - this documentation, please review our style guidelines - (below). When editing the source, please refrain from using - elements with the wrong semantic meaning for aesthetic - reasons. Aesthetic issues can be addressed by changes to the theme. +.. code-block:: python -For **titles and headers**: + import gcpy + print("hello world") -* Section headers should be underlined by :literal:`#` characters -* Subsection headers should be underlined by :literal:`-` characters -* Subsubsection headers should be underlined by :literal:`^` characters -* Subsubsubsection headers should be avoided, but if necessary, they should be underlined by :literal:`"` characters +Fortran +~~~~~~~ + +.. code-block:: none -**File paths** (including directories) occuring in the text should use -the :literal:`:file:` role. + .. code-block:: Fortran -**Program names** (e.g. :program:`cmake`) occuring in the text should -use the :literal:`:program:` role. + DO I = 1, 10 + PRINT*, I + ENDDO -**OS-level commands** (e.g. :command:`rm`) occuring in the text should -use the :literal:`:command:` role. +Renders as: -**Environment variables** occuring in the text should use the -:literal:`:envvar:` role. +.. code-block:: Fortran -**Inline code** or code variables occuring in the text should use the -:literal:`:code:` role. + DO I = 1, 10 + PRINT*, I + ENDDO -**Code snippets** should use :literal:`.. code-block:: ` -directive like so +Bash +~~~~ .. code-block:: none - .. code-block:: python + .. code-block:: bash - import gcpy - print("hello world") + #!/bin/bash + + for f in *.nc; do + echo $f + done + +Renders as: + +.. code-block:: bash + + #!/bin/bash + + for f in *.nc; do + echo $f + done -The language can be "none" to omit syntax highlighting. +Command line (aka "console") +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -For command line instructions, the "console" language should be -used. The :literal:`$` should be used to denote the console's -prompt. If the current working directory is relevant to the -instructions, a prompt like :literal:`gcuser:~/path1/path2$` should be -used. +.. code-block:: none + + .. code-block:: console + + $ ls -l $HOME + +Renders as: + +.. code-block:: console + + $ ls -l $HOME + +No formatting +~~~~~~~~~~~~~ + +.. code-block:: none + + .. code-block:: none + + This text renders without any special formatting. + +Renders as: + +.. code-block:: none -**Inline literals** (e.g. the :literal:`$` above) should use the -:literal:`:literal:` role. + This text renders without any special formatting. diff --git a/docs/source/tech_info/06_info_for_kpp_developers.rst b/docs/source/tech_info/06_info_for_kpp_developers.rst index ab21bc3..8c513dc 100644 --- a/docs/source/tech_info/06_info_for_kpp_developers.rst +++ b/docs/source/tech_info/06_info_for_kpp_developers.rst @@ -18,9 +18,10 @@ KPP directory structure The KPP distribution will unfold a directory :envvar:`$KPP_HOME` with the following subdirectories: -.. option:: src/ +src/ +---- - Contains the KPP source code files: +Contains the KPP source code files: .. _table-kpp-dirs: @@ -67,72 +68,87 @@ following subdirectories: | :file:`y.tab.h` | Bison generated header file | +-----------------------+-------------------------------------+ -.. option:: bin/ +bin/ +---- - Contains the KPP executable. This directory should be added to the - :envvar:`PATH` environment variable. +Contains the KPP executable. This directory should be added to the +:envvar:`PATH` environment variable. -.. option:: util/ +util/ +----- - Contains different function templates useful for the simulation. Each - template file has a suffix that matches the appropriate target - language (Fortran90, C, or Matlab). KPP will run the template files - through the substitution preprocessor (cf. - :ref:`list-of-symbols-replaced`). The user can define their own - auxiliary functions by inserting them into the files. +Contains different function templates useful for the simulation. Each +template file has a suffix that matches the appropriate target +language (Fortran90, C, or Matlab). KPP will run the template files +through the substitution preprocessor (cf. +:ref:`list-of-symbols-replaced`). The user can define their own +auxiliary functions by inserting them into the files. -.. option:: models/ +models/ +------- - Contains the description of the chemical models. Users - can define their own models by placing the model description files in - this directory. The KPP distribution contains several models from - atmospheric chemistry which can be used as templates for model - definitions. +Contains the description of the chemical models. Users +can define their own models by placing the model description files in +this directory. The KPP distribution contains several models from +atmospheric chemistry which can be used as templates for model +definitions. -.. option:: drv/ +drv/ +---- - Contains driver templates for chemical simulations. Each driver has a - suffix that matches the appropriate target language (Fortran90, C, or - Matlab). KPP will run the appropriate driver through the substitution - preprocessor (cf. :ref:`list-of-symbols-replaced`). Users can also - define their own driver templates here. +Contains driver templates for chemical simulations. Each driver has a +suffix that matches the appropriate target language (Fortran90, C, or +Matlab). KPP will run the appropriate driver through the substitution +preprocessor (cf. :ref:`list-of-symbols-replaced`). Users can also +define their own driver templates here. -.. option:: int/ +int/ +---- - Contains numerical solvers (integrators). The :command:`#INTEGRATOR` - command will force KPP to look into this directory for a definition - file with suffix :code:`.def`. This file selects the numerical solver - etc. Each integrator template is found in a file that ends with the - appropriate suffix (:code:`.f90`, :code:`.c`, or :code:`.m`). The - selected template is processed by the substitution preprocessor (cf. - :ref:`list-of-symbols-replaced`). Users can define their own - numerical integration routines in the :code:`user_contributed` - subdirectory. +Contains numerical solvers (integrators). The :command:`#INTEGRATOR` +command will force KPP to look into this directory for a definition +file with suffix :code:`.def`. This file selects the numerical solver +etc. Each integrator template is found in a file that ends with the +appropriate suffix (:code:`.f90`, :code:`.c`, or :code:`.m`). The +selected template is processed by the substitution preprocessor (cf. +:ref:`list-of-symbols-replaced`). Users can define their own +numerical integration routines in the :code:`user_contributed` +subdirectory. -.. option:: examples/ +examples/ +--------- - Contains several model description examples (:file:`.kpp` files) - which can be used as templates for building simulations with KPP. +Contains several model description examples (:file:`.kpp` files) +which can be used as templates for building simulations with KPP. -.. option:: site-lisp/ +site-lisp/ +---------- +Contains the file :file:`kpp.el` which provides a KPP mode for emacs +with color highlighting. - Contains the file :file:`kpp.el` which provides a KPP mode for emacs - with color highlighting. +ci-tests/ +--------- -.. option:: ci-tests/ +Contains directories defining several :ref:`ci-tests`. - Contains directories defining several :ref:`ci-tests`. +.ci-pipelines/ +-------------- -.. option:: .ci-pipelines/ +Hidden directory containing a YAML file with settings for automatically +running the continuous integration tests on `Azure DevOps Pipelines +`_ - Hidden directory containing a YAML file with settings for automatically - running the continuous integration tests on `Azure DevOps Pipelines - `_ +Also contains bash scripts (ending in :file:`.sh`) for running the +continuous integration tests either automatically in Azure Dev +Pipelines, or manually from the command line. For more +information, please see :ref:`ci-tests`. - Also contains bash scripts (ending in :file:`.sh`) for running the - continuous integration tests either automatically in Azure Dev - Pipelines, or manually from the command line. For more - information, please see :ref:`ci-tests`. +.github/workflows +----------------- + +Contains configuration files for `GitHub Actions +`_ that will run +automatically when commits are pushed or when pull requests are opened. .. _kpp-env-vars: @@ -151,42 +167,47 @@ The :file:`$KPP_HOME/bin` directory. should be added to the There are also several optional environment variable that control the places where KPP looks for module files, integrators, and drivers: -.. option:: KPP_HOME +KPP_HOME +-------- - Required, stores the absolute path to the KPP distribution. +Required, stores the absolute path to the KPP distribution. - Default setting: none. +Default setting: none. -.. option:: KPP_FLEX_LIB_DIR +KPP_FLEX_LIB_DIR +---------------- - Optional. Use this to specify the path to the :ref:`flex library - file ` (:file:`libfl.so` or :file:`libfl.a`) that are - needed to :ref:`build the KPP executable `. The KPP - build sequence will use the path contained in - :envvar:`KPP_FLEX_LIB_DIR` if the flex library file cannot be found - in :file:`/usr/lib`, :file:`/usr/lib64`, and similar standard - library paths. +Optional. Use this to specify the path to the :ref:`flex library +file ` (:file:`libfl.so` or :file:`libfl.a`) that are +needed to :ref:`build the KPP executable `. The KPP +build sequence will use the path contained in +:envvar:`KPP_FLEX_LIB_DIR` if the flex library file cannot be found +in :file:`/usr/lib`, :file:`/usr/lib64`, and similar standard +library paths. -.. option:: KPP_MODEL +KPP_MODEL +--------- - Optional, specifies additional places where KPP will look for model - files before searching the default location. +Optional, specifies additional places where KPP will look for model +files before searching the default location. - Default setting: :file:`$KPP_HOME/models`. +Default setting: :file:`$KPP_HOME/models`. -.. option:: KPP_INT +KPP_INT +------- - Optional, specifies additional places where KPP will look for - integrator files before searching the default. +Optional, specifies additional places where KPP will look for +integrator files before searching the default. - Default setting: :file:`$KPP_HOME/int`. +Default setting: :file:`$KPP_HOME/int`. -.. option:: KPP_DRV +KPP_DRV +------- - Optional specifies additional places where KPP will look for driver - files before searching the default directory. +Optional specifies additional places where KPP will look for driver +files before searching the default directory. - Default setting: :file:`$KPP_HOME/drv`. +Default setting: :file:`$KPP_HOME/drv`. .. _kpp-internal-modules: @@ -330,13 +351,12 @@ then runs a short "box model" simulation with the generated code. C-I tests help to ensure that new features and updates added to KPP will not break any existing functionality. -The continuous integration tests will run automatically on `Azure -DevOps Pipelines -`_ each time a -commit is pushed to the `KPP Github repository -`_. You can also run the -integration tests :ref:`locally on your own computer -`. +C-I tests will run automatically as a `GitHub Action +`_ when commits +are pushed to the `KPP Github repository +`_, or when a new pull +requests are opened. You may also run the integration tests +:ref:`locally on your own computer `. .. _list-of-ci-tests: @@ -395,6 +415,10 @@ List of continuous integration tests - Fortran90 - small_strato - runge_kutta + * - :code:`F90_rkadj` + - Fortran90 + - small_strato + - runge_kutta_adj * - :code:`F90_rktlm` - Fortran90 - small_strato @@ -435,6 +459,10 @@ List of continuous integration tests - Fortran90 - saprcnov - rosenbrock + * - :code:`F90_sd4` + - Fortran90 + - small_strato + - sdirk4 * - :code:`F90_sd` - Fortran90 - small_strato @@ -443,6 +471,10 @@ List of continuous integration tests - Fortran90 - small_strato - sdirk_adj + * - :code:`F90_sdtlm` + - Fortran90 + - small_strato + - sdirk_tlm * - :code:`F90_seulex` - Fortran90 - saprcnov @@ -464,84 +496,64 @@ Notes about C-I tests: :ref:`running-kpp-with-an-example-mechanism`. #. :file:`X_minver` tests if the :ref:`minversion-cmd` command works properly. -#. Due to memory restrictions, the :file:`F90_mcm` and - :file:`F90_mcm_h211b` are not run on the Microsoft Azure Dev - pipelines platform. However, you can run these tests manaully. Each continuous integration test is contained in a subdirectory of :file:`$KPP_HOME/ci-tests`. In each subdirectory is a KPP definition file (ending in :file:`.kpp`). -.. _running-ci-tests-on-azure: +.. _running-ci-tests: -Running continuous integration tests on Azure DevOps Pipelines --------------------------------------------------------------- +Running continuous integration tests as a GitHub Action +------------------------------------------------------- -The files needed to run the C-I tests are located in the -:file:`$KPP_HOME/.ci-pipelines` directory: +The files needed to run the C-I tests are described below. -.. _table-ci-pipelines: +.. _running-ci-tests-action: -.. table:: Files needed to execute C-I tests - :align: center +run-ci-tests.yml +~~~~~~~~~~~~~~~~ + +**Path**: :file:`$KPP_HOME/.github/workflows/run-ci-tests.yml` + +**Description:** Configuration file with commands to download KPP, load +libraries, and run the C-I tests as a GitHub Action. + +C-I tests will run automatically when a commit is pushed to any branch +at `https://github.com/KineticPreProcessor/KPP +`_, or when a new pull +request is opened there. This is the recommended setting, but you can +restrict this so that only pushes or pull requests to certain branches +will trigger the C-I tests. + +.. _running-ci-tests-testing: + +ci-testing-script.sh +~~~~~~~~~~~~~~~~~~~~ + +**Path:** :file:`$KPP_HOME/.ci-pipelines/ci-testing-script.sh` + +**Description:** Runs the KPP C-I tests as a GitHub Action, or on a +local computer system. + +.. _running-ci-tests-cleanup: + +ci-cleanup-script.sh +~~~~~~~~~~~~~~~~~~~~ + +**Path:** :file:`$KPP_HOME/.ci-pipelines/ci-cleanup-script.sh` + +**Description:** Removes compiler-generated files (e.g. :file:`*.o`, +:file:`.mod` , and :file:`.exe`) from C-I test folders. + +.. _running-ci-tests-defs: + +ci-common-defs.sh +~~~~~~~~~~~~~~~~~~ + +**Path:** :file:`$KPP_HOME/.ci-pipelines/ci-common-defs.sh` - +-------------------------------+------------------------------------------+ - | File | Description | - +===============================+==========================================+ - | :file:`Dockerfile` | File containing specifications for the | - | | Docker container that will be used to | - | | run C-I tests on Azure DevOps Pipelines. | - | | Also contains commands needed to run | - | | the C-I scripts in the Docker container. | - +-------------------------------+------------------------------------------+ - | :file:`build_testing.yml` | Contains options for triggering C-I | - | | tests on Azure DevOps Pipelines. | - +-------------------------------+------------------------------------------+ - | :file:`ci-testing-script.sh` | Driver script for running C-I tests. | - | | Can be used on Azure DevOps Pipelines | - | | or on a local computer. | - +-------------------------------+------------------------------------------+ - | :file:`ci-cleanup-script.sh` | Script to remove compiler-generated | - | | files (e.g. ``*.o``, ``.mod``, and | - | | ``.exe``) from C-I test folders. | - +-------------------------------+------------------------------------------+ - | :file:`ci-common-defs.sh` | Script with common variable and function | - | | definitions needed by | - | | :file:`ci-testing-script.sh` and | - | | :file:`ci-cleanup-script.sh`. | - +-------------------------------+------------------------------------------+ - -The :file:`Dockerfile` contains the software environment for `Azure -DevOps Pipelines -`_. You -should not have to update this file. - -File :file:`build_testing.yml` defines the runtime options for Azure -DevOps Pipelines. The following settings determine which branches -will trigger C-I tests: - -.. code-block:: yaml - - # Run a C-I test when a push to any branch is made. - trigger: - branches: - include: - - '*' - pr: - branches: - include: - - '*' - -Currently this is set to trigger the C-I tests when a commit or pull -request is made to any branch of -`https://github.com/KineticPreProcessor/KPP -`_. This is the recommended -setting, but you can restrict this so that only pushes or pull requests -to certain branches will trigger the C-I tests. - -The script :file:`ci-testing-script.sh` executes all of the C-I tests -whenever a push or a pull request is made to the selected branches in -the KPP Github repository. +**Description** Contains common variable and function definitions needed by +:ref:`running-ci-tests-testing` and :ref:`running-ci-tests-cleanup`. .. _running-ci-tests-locally: diff --git a/docs/source/tech_info/07_numerical_methods.rst b/docs/source/tech_info/07_numerical_methods.rst index 66c40a2..88b0a27 100644 --- a/docs/source/tech_info/07_numerical_methods.rst +++ b/docs/source/tech_info/07_numerical_methods.rst @@ -380,7 +380,7 @@ Rosenbrock with H211b time stepping **Integrator file:** :file:`int/rosenbrock_h211b_qssa.f90` -H211b time stepping according to :cite:t:`Soederlind_2022`, as +H211b time stepping according to :cite:t:`Soederlind_2003`, as implemented by :cite:t:`Dreger_2025`. .. _rk-methods: diff --git a/examples/rkadj.kpp b/examples/rkadj.kpp new file mode 100644 index 0000000..3e28e80 --- /dev/null +++ b/examples/rkadj.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR runge_kutta_adj +#LANGUAGE Fortran90 +#DRIVER general_adj diff --git a/examples/sd4.kpp b/examples/sd4.kpp new file mode 100644 index 0000000..7dc9ea9 --- /dev/null +++ b/examples/sd4.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR sdirk4 +#LANGUAGE Fortran90 +#DRIVER general diff --git a/examples/sdtlm.kpp b/examples/sdtlm.kpp new file mode 100644 index 0000000..9865d6b --- /dev/null +++ b/examples/sdtlm.kpp @@ -0,0 +1,4 @@ +#MODEL small_strato +#INTEGRATOR sdirk_tlm +#LANGUAGE Fortran90 +#DRIVER general_tlm diff --git a/int/dvode.f90 b/int/dvode.f90 index 3e5f7ec..64430e4 100644 --- a/int/dvode.f90 +++ b/int/dvode.f90 @@ -4,7 +4,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP - USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, Set2zero, WLAMCH + USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve IMPLICIT NONE PUBLIC diff --git a/int/lsode.f90 b/int/lsode.f90 index 93a4bb8..3b6d73f 100644 --- a/int/lsode.f90 +++ b/int/lsode.f90 @@ -11,7 +11,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, WLAMCH + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC diff --git a/int/radau5.f90 b/int/radau5.f90 index 89a4125..f409dde 100644 --- a/int/radau5.f90 +++ b/int/radau5.f90 @@ -412,7 +412,7 @@ SUBROUTINE RADAU5(N,T,Tend,Y,H,RelTol,AbsTol, & !~~~> Roundoff SMALLEST NUMBER SATISFYING 1.0d0+Roundoff>1.0d0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> RCNTRL(1) = Hmin - not used Hmin = ZERO @@ -642,12 +642,12 @@ SUBROUTINE RAD_Integrator( & ! STARTING VALUES FOR NEWTON ITERATION !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) - CALL Set2zero(N,F1) - CALL Set2zero(N,F2) - CALL Set2zero(N,F3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp + F1(1:N) = 0.0_dp + F2(1:N) = 0.0_dp + F3(1:N) = 0.0_dp ELSE C3Q=H/Hold C1Q=rkC(1)*C3Q @@ -726,9 +726,9 @@ SUBROUTINE RAD_Integrator( & END IF END IF DYNOLD=MAX(DYNO,Roundoff) - CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1 - CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2 - CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3 + F1(1:N) = F1(1:N) + Z1(1:N) ! F1 <- F1 + Z1 + F2(1:N) = F2(1:N) + Z2(1:N) ! F2 <- F2 + Z2 + F3(1:N) = F3(1:N) + Z3(1:N) ! F3 <- F3 + Z3 ! Z(1,2,3) = Transf x F(1,2,3) CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3) NewtonDone = (FacConv*DYNO <= TolNewton) diff --git a/int/rosenbrock.f90 b/int/rosenbrock.f90 index e5ae740..2460c8d 100644 --- a/int/rosenbrock.f90 +++ b/int/rosenbrock.f90 @@ -335,7 +335,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -519,9 +519,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -581,55 +578,49 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & RETURN END IF -!~~~> Compute the stages -Stage: DO istage = 1, ros_S + !~~~> Compute the stages + Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N) - ioffset = N*(istage-1) + ioffset = N*(istage-1) ! For the 1st istage the function has been computed previously - IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) + IF ( istage == 1 ) THEN Fcn(1:N) = Fcn0(1:N) - ! istage>1 and a new function evaluation is needed at the current istage - ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) + ! istage>1 and a new function evaluation is needed at the current istage + ELSEIF ( ros_NewF(istage) ) THEN Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) & + + ros_A((istage-1)*(istage-2)/2+j) * K(N*(j-1)+1:N*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) - K(ioffset+1:ioffset+N) = Fcn(1:N) - DO j = 1, istage-1 + END IF ! if istage == 1 elseif ros_NewF(istage) + K(ioffset+1:ioffset+N) = Fcn(1:N) + DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) - END DO - IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HC * K(N*(j-1)+1:N*j) + END DO + IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) - END IF - CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HG * dFdT(1:N) + END IF + CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) END DO Stage - !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -645,7 +636,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -732,8 +722,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -781,16 +771,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_adj.f90 b/int/rosenbrock_adj.f90 index 90c6e2b..aa644f9 100644 --- a/int/rosenbrock_adj.f90 +++ b/int/rosenbrock_adj.f90 @@ -438,7 +438,7 @@ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, & !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -747,8 +747,6 @@ SUBROUTINE ros_DPush( S, T, H, Ystage, K, E, P ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - !CALL WCOPY(NVAR*S,Ystage,1,chk_Y(1,stack_ptr),1) - !CALL WCOPY(NVAR*S,K,1,chk_K(1,stack_ptr),1) chk_Y(1:NVAR*S,stack_ptr) = Ystage(1:NVAR*S) chk_K(1:NVAR*S,stack_ptr) = K(1:NVAR*S) IF (SaveLU) THEN @@ -784,11 +782,8 @@ SUBROUTINE ros_DPop( S, T, H, Ystage, K, E, P ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - !CALL WCOPY(NVAR*S,chk_Y(1,stack_ptr),1,Ystage,1) - !CALL WCOPY(NVAR*S,chk_K(1,stack_ptr),1,K,1) Ystage(1:NVAR*S) = chk_Y(1:NVAR*S,stack_ptr) K(1:NVAR*S) = chk_K(1:NVAR*S,stack_ptr) - !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1) IF (SaveLU) THEN #ifdef FULL_ALGEBRA E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr) @@ -817,9 +812,6 @@ SUBROUTINE ros_CPush( T, H, Y, dY, d2Y ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - !CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - !CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1) - !CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1) chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR) chk_d2Y(1:NVAR,stack_ptr) = d2Y(1:NVAR) @@ -840,9 +832,6 @@ SUBROUTINE ros_CPop( T, H, Y, dY, d2Y ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - !CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) - !CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1) - !CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr) d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr) @@ -1004,60 +993,59 @@ SUBROUTINE ros_FwdInt ( Y, & RETURN END IF -!~~~> Compute the stages -Stage: DO istage = 1, ros_S + !~~~> Compute the stages + Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR) - ioffset = NVAR*(istage-1) + ioffset = NVAR*(istage-1) ! For the 1st istage the function has been computed previously - IF ( istage == 1 ) THEN - CALL WCOPY(NVAR,Fcn0,1,Fcn,1) + IF ( istage == 1 ) THEN + Fcn(1:NVAR) = Fcn0(1:NVAR) IF (AdjointType == Adj_discrete) THEN ! Save stage solution - ! CALL WCOPY(NVAR,Y,1,Ystage(1),1) Ystage(1:NVAR) = Y(1:NVAR) - CALL WCOPY(NVAR,Y,1,Ynew,1) + Ynew(1:NVAR) = Y(1:NVAR) END IF - ! istage>1 and a new function evaluation is needed at the current istage - ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR,Y,1,Ynew,1) + ! istage>1 and a new function evaluation is needed at the current istage + ELSEIF ( ros_NewF(istage) ) THEN + Ynew(1:NVAR) = Y(1:NVAR) DO j = 1, istage-1 - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1),1,Ynew,1) + Ynew(1:NVAR) = Ynew(1:NVAR) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - END IF ! if istage == 1 elseif ros_NewF(istage) + END IF ! if istage == 1 elseif ros_NewF(istage) ! save stage solution every time even if ynew is not updated - IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN - ! CALL WCOPY(NVAR,Ynew,1,Ystage(ioffset+1),1) + IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR) - END IF - CALL WCOPY(NVAR,Fcn,1,K(ioffset+1),1) - DO j = 1, istage-1 + END IF + K(ioffset+1:ioffset+NVAR) = Fcn(1:NVAR) + DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1) - END DO - IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN + K(ioffset+1:ioffset+NVAR) = K(ioffset+1:ioffset+NVAR) & + + HC * K(NVAR*(j-1)+1:NVAR*j) + END DO + IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1) - END IF - CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1)) + K(ioffset+1:ioffset+NVAR) = K(ioffset+1:ioffset+NVAR) + HG * dFdT(1:NVAR) + END IF + CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1)) END DO Stage - !~~~> Compute the new solution - CALL WCOPY(NVAR,Y,1,Ynew,1) + Ynew(1:NVAR) = Y(1:NVAR) DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1) + Ynew(1:NVAR) = Ynew(1:NVAR) + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j) END DO !~~~> Compute the error estimation - CALL WSCAL(NVAR,ZERO,Yerr,1) + Yerr(1:NVAR) = ZERO DO j=1,ros_S - CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1) + Yerr(1:NVAR) = Yerr(1:NVAR) + ros_E(j) * K(NVAR*(j-1)+1:NVAR*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -1079,11 +1067,11 @@ SUBROUTINE ros_FwdInt ( Y, & CALL Jac_SP_Vec( Jac0, Fcn0, K(1) ) #endif IF (.NOT. Autonomous) THEN - CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) + K(1:NVAR) = K(1:NVAR) + dFdT(1:NVAR) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) END IF - CALL WCOPY(NVAR,Ynew,1,Y,1) + Y(1:NVAR) = Ynew(1:NVAR) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1126,7 +1114,7 @@ SUBROUTINE ros_FwdInt ( Y, & #endif IF (.NOT. Autonomous) THEN CALL ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) - CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) + K(1:NVAR) = K(1:NVAR) + dFdT(1:NVAR) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) !~~~> Deallocate stage buffer: only needed for discrete adjoint @@ -1212,7 +1200,7 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & Ghimj(i,i) = Ghimj(i,i)+Tau END DO #else - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau END DO @@ -1231,16 +1219,18 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & !~~~> Compute U DO m = 1,NADJ - CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1) - CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1) + U(istart:istart+NVAR-1,m) = Lambda(1:NVAR,m) + U(istart:istart+NVAR-1,m) = ros_M(istage) * U(istart:istart+NVAR-1,m) END DO ! m=1:NADJ DO j = istage+1, ros_S jstart = NVAR*(j-1) + 1 HA = ros_A((j-1)*(j-2)/2+istage) HC = ros_C((j-1)*(j-2)/2+istage)/(Direction*H) DO m = 1,NADJ - CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1) - CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1) + U(istart:istart+NVAR-1,m) = U(istart:istart+NVAR-1,m) & + + HA * V(jstart:jstart+NVAR-1,m) + U(istart:istart+NVAR-1,m) = U(istart:istart+NVAR-1,m) & + + HC * U(jstart:jstart+NVAR-1,m) END DO ! m=1:NADJ END DO DO m = 1,NADJ @@ -1273,10 +1263,10 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & istart = NVAR*(istage-1) + 1 DO m = 1,NADJ ! Add V_i - CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + V(istart:istart+NVAR-1,m) ! Add (H0xK_i)^T * U_i CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp ) - CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + Tmp(1:NVAR) END DO ! m=1:NADJ END DO ! Add H * dJac_dT_0^T * \sum(gamma_i U_i) @@ -1286,14 +1276,15 @@ SUBROUTINE ros_DadjInt ( NADJ, Lambda, & Tmp(1:NVAR) = ZERO DO istage = 1, ros_S istart = NVAR*(istage-1) + 1 - CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1) + Tmp(1:NVAR) = Tmp(1:NVAR) & + + ros_Gamma(istage) * U(istart:istart+NVAR-1,m) END DO #ifdef FULL_ALGEBRA Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp) #else CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2) #endif - CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1) + Lambda(1:NVAR,m) = Lambda(1:NVAR,m) + H * Tmp2(1:NVAR) END DO ! m=1:NADJ END IF ! .NOT.Autonomous @@ -1404,7 +1395,7 @@ SUBROUTINE ros_CadjInt ( & #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif - CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) + dFdT(1:NVAR,iadj) = -dFdT(1:NVAR,iadj) END DO END IF @@ -1412,7 +1403,7 @@ SUBROUTINE ros_CadjInt ( & #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) + Jac0(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1441,16 +1432,17 @@ SUBROUTINE ros_CadjInt ( & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) + Fcn(1:NVAR,iadj) = Fcn0(1:NVAR,iadj) END DO ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) + Ynew(1:NVAR,1:NADJ) = Y(1:NVAR,1:NADJ) DO j = 1, istage-1 - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) - END DO + DO iadj = 1, NADJ + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO Tau = T + ros_Alpha(istage)*Direction*H ! CALL FunTemplate(Tau,Ynew,Fcn) @@ -1461,7 +1453,7 @@ SUBROUTINE ros_CadjInt ( & #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) + Jac(1:LU_NONZERO) = -Jac(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1469,25 +1461,25 @@ SUBROUTINE ros_CadjInt ( & #else CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj)) #endif - !CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1) END DO END IF ! if istage == 1 elseif ros_NewF(istage) DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) + K(ioffset+1:ioffset+NVAR,iadj) = Fcn(1:NVAR,iadj) END DO DO j = 1, istage-1 - HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & - K(ioffset+1,iadj),1) - END DO + HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HC * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN - HG = Direction*H*ros_Gamma(istage) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) - END DO + HG = Direction*H*ros_Gamma(istage) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HG * dFdT(1:NVAR,iadj) + END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) @@ -1498,18 +1490,20 @@ SUBROUTINE ros_CadjInt ( & !~~~> Compute the new solution DO iadj = 1, NADJ - CALL WCOPY(NVAR,Y(1,iadj),1,Ynew(1,iadj),1) + Ynew(1:NVAR,iadj) = Y(1:NVAR,iadj) DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) END DO END DO !~~~> Compute the error estimation - CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1) + Yerr(1:NVAR,1:NADJ) = ZERO DO j=1,ros_S - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1,iadj),1,Yerr(1,iadj),1) - END DO + DO iadj = 1, NADJ + Yerr(1:NVAR,iadj) = Yerr(1:NVAR,iadj) & + + ros_E(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO !~~~> Max error among all adjoint components iadj = 1 @@ -1524,7 +1518,7 @@ SUBROUTINE ros_CadjInt ( & ! ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1) + Y(1:NVAR,1:NADJ) = Ynew(1:NVAR,1:NADJ) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1598,9 +1592,6 @@ SUBROUTINE ros_SimpleCadjInt ( & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1617,7 +1608,6 @@ SUBROUTINE ros_SimpleCadjInt ( & T = chk_T(istack) H = chk_H(istack-1) - !CALL WCOPY(NVAR,chk_Y(1,istack),1,Y0,1) Y0(1:NVAR) = chk_Y(1:NVAR,istack) !~~~> Compute the Jacobian at current time @@ -1633,7 +1623,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif - CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) + dFdT(1:NVAR,iadj) = -dFdT(1:NVAR,iadj) END DO END IF @@ -1641,7 +1631,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) + Jac0(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1659,8 +1649,7 @@ SUBROUTINE ros_SimpleCadjInt ( & Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv END DO @@ -1681,17 +1670,16 @@ SUBROUTINE ros_SimpleCadjInt ( & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) - END DO + Fcn(1:NVAR,1:NADJ) = Fcn0(1:NVAR,1:NADJ) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) + Ynew(1:NVAR,1:NADJ) = Y(1:NVAR,1:NADJ) DO j = 1, istage-1 - DO iadj = 1, NADJ - CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & - K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) - END DO + DO iadj = 1, NADJ + Ynew(1:NVAR,iadj) = Ynew(1:NVAR,iadj) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO Tau = T + ros_Alpha(istage)*Direction*H CALL ros_Hermite3( chk_T(istack-1), chk_T(istack), Tau, & @@ -1702,7 +1690,7 @@ SUBROUTINE ros_SimpleCadjInt ( & #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else - CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) + Jac(1:LU_NONZERO) = -Jac(1:LU_NONZERO) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA @@ -1712,22 +1700,20 @@ SUBROUTINE ros_SimpleCadjInt ( & #endif END DO END IF ! if istage == 1 elseif ros_NewF(istage) - - DO iadj = 1, NADJ - CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) - END DO + K(ioffset+1:ioffset+NVAR,1:NADJ) = Fcn(1:NVAR,1:NADJ) DO j = 1, istage-1 - HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & - K(ioffset+1,iadj),1) - END DO + HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HC * K(NVAR*(j-1)+1:NVAR*j,iadj) + END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN - HG = Direction*H*ros_Gamma(istage) - DO iadj = 1, NADJ - CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) - END DO + HG = Direction*H*ros_Gamma(istage) + DO iadj = 1, NADJ + K(ioffset+1:ioffset+NVAR,iadj) = K(ioffset+1:ioffset+NVAR,iadj) & + + HG * dFdT(1:NVAR,iadj) + END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) @@ -1739,7 +1725,8 @@ SUBROUTINE ros_SimpleCadjInt ( & !~~~> Compute the new solution DO iadj = 1, NADJ DO j=1,ros_S - CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Y(1,iadj),1) + Y(1:NVAR,iadj) = Y(1:NVAR,iadj) & + + ros_M(j) * K(NVAR*(j-1)+1:NVAR*j,iadj) END DO END DO @@ -1801,8 +1788,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(NVAR,(ONE/Delta),dFdT,1) + dFdT(1:NVAR) = dFdT(1:NVAR) - Fcn0(1:NVAR) + dFdT(1:NVAR) = dFdT(1:NVAR) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1830,11 +1817,11 @@ SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, Jac0, dJdT ) CALL JacTemplate( T+Delta, Y, dJdT ) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA - CALL WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1) + dJdT(1:NVAR,1:NVAR) = dJdT(1:NVAR,1:NVAR) - Jac0(1:NVAR,1:NVAR) + dJdT(1:NVAR,1:NVAR) = dJdT(1:NVAR,1:NVAR) * (ONE/Delta) #else - CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) - Jac0(1:LU_NONZERO) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) * (ONE/Delta) #endif END SUBROUTINE ros_JacTimeDerivative @@ -1884,15 +1871,15 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - CALL WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1) - CALL WSCAL(NVAR*NVAR,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv @@ -2058,30 +2045,30 @@ SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y ) amb(i) = amb(i-1)*amb(1) END DO - ! c(1) = ya; - CALL WCOPY(NVAR,Ya,1,C(1,1),1) + C(1:NVAR,1) = Ya(1:NVAR) + ! c(2) = ja; - CALL WCOPY(NVAR,Ja,1,C(1,2),1) + C(1:NVAR,2) = Ja(1:NVAR) + ! c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb ; - CALL WCOPY(NVAR,Ya,1,C(1,3),1) - CALL WSCAL(NVAR,-3.0*amb(2),C(1,3),1) - CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1,3),1) - CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1,3),1) - CALL WAXPY(NVAR,amb(1),Jb,1,C(1,3),1) + C(1:NVAR,3) = 2.0*amb(1) * Ja(1:NVAR) & + + amb(1) * Jb(1:NVAR) & + - 3.0*amb(2) * Ya(1:NVAR) & + + 3.0*amb(2) * Yb(1:NVAR) + ! c(4) = 1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb ; - CALL WCOPY(NVAR,Ya,1,C(1,4),1) - CALL WSCAL(NVAR,-2.0*amb(3),C(1,4),1) - CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1,4),1) - CALL WAXPY(NVAR,amb(2),Ja,1,C(1,4),1) - CALL WAXPY(NVAR,amb(2),Jb,1,C(1,4),1) + C(1:NVAR,4) = amb(2) * Ja(1:NVAR) & + + amb(2) * Jb(1:NVAR) & + - 2.0*amb(3) * Ya(1:NVAR) & + + 2.0*amb(3) * Yb(1:NVAR) +! Unrolled loop: Y = Tau^3*c(4) + Tau^2*c(3) + Tau*c(2) + c(1) Tau = T - a - CALL WCOPY(NVAR,C(1,4),1,Y,1) - CALL WSCAL(NVAR,Tau**3,Y,1) - DO j = 3,1,-1 - CALL WAXPY(NVAR,TAU**(j-1),C(1,j),1,Y,1) - END DO + Y(1:NVAR) = Tau**3 * C(1:NVAR,4) & + + Tau**2 * C(1:NVAR,3) & + + Tau * C(1:NVAR,2) & + + C(1:NVAR,1) END SUBROUTINE ros_Hermite3 @@ -2109,45 +2096,43 @@ SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y ) END DO ! c(1) = ya; - CALL WCOPY(NVAR,Ya,1,C(1,1),1) + C(1:NVAR,1) = Ya(1:NVAR) + ! c(2) = ja; - CALL WCOPY(NVAR,Ja,1,C(1,2),1) + C(1:NVAR,2) = Ja(1:NVAR) + ! c(3) = ha/2; - CALL WCOPY(NVAR,Ha,1,C(1,3),1) - CALL WSCAL(NVAR,HALF,C(1,3),1) + C(1:NVAR,3) = Ha(1:NVAR) * HALF ! c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb + 1.5*amb(1)*ha - 0.5*amb(1)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,4),1) - CALL WSCAL(NVAR,10.0*amb(3),C(1,4),1) - CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1,4),1) - CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1,4),1) - CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1,4),1) - CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1,4),1) - CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1,4),1) + C(1:NVAR,4) = 10.0*amb(3) * Ya(1:NVAR) & + - 10.0*amb(3) * Yb(1:NVAR) & + - 6.0*amb(2) * Ja(1:NVAR) & + - 4.0*amb(2) * Jb(1:NVAR) & + + 1.5*amb(1) * Ha(1:NVAR) & + - 0.5*amb(1) * Hb(1:NVAR) ! c(5) = 15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb + 1.5*amb(2)*ha - 1*amb(2)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,5),1) - CALL WSCAL(NVAR, 15.0*amb(4),C(1,5),1) - CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1,5),1) - CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1,5),1) - CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1,5),1) - CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1,5),1) - CALL WAXPY(NVAR,-amb(2),Hb,1,C(1,5),1) + C(1:NVAR,5) = 15.0*amb(4) * Ya(1:NVAR) & + - 15.0*amb(4) * Yb(1:NVAR) & + - 8.0*amb(3) * Ja(1:NVAR) & + - 7.0*amb(3) * Jb(1:NVAR) & + + 1.5*amb(2) * Ha(1:NVAR) & + - amb(2) * Hb(1:NVAR) ! c(6) = 6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb + 0.5*amb(3)*ha -0.5*amb(3)*hb ; - CALL WCOPY(NVAR,Ya,1,C(1,6),1) - CALL WSCAL(NVAR, 6.0*amb(5),C(1,6),1) - CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1,6),1) - CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1,6),1) - CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1,6),1) - CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1,6),1) - CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1,6),1) + C(1:NVAR,6) = 6.0*amb(5) * Ya(1:NVAR) & + - 6.0*amb(5) * Yb(1:NVAR) & + - 3.0*amb(4) * Ja(1:NVAR) & + - 3.0*amb(4) * Jb(1:NVAR) & + + 0.5*amb(3) * Ha(1:NVAR) & + - 0.5*amb(3) * Hb(1:NVAR) + Tau = T - a - CALL WCOPY(NVAR,C(1,6),1,Y,1) + Y(1:NVAR) = C(1:NVAR,6) DO j = 5,1,-1 - CALL WSCAL(NVAR,Tau,Y,1) - CALL WAXPY(NVAR,ONE,C(1,j),1,Y,1) + Y(1:NVAR) = Tau * Y(1:NVAR) + C(1:NVAR,j) END DO END SUBROUTINE ros_Hermite5 diff --git a/int/rosenbrock_autoreduce.f90 b/int/rosenbrock_autoreduce.f90 index a7a8e07..b6e2821 100644 --- a/int/rosenbrock_autoreduce.f90 +++ b/int/rosenbrock_autoreduce.f90 @@ -307,7 +307,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ISTATUS(1:8) = 0 RSTATUS(1:4) = ZERO -!~~~> Autonomous or time dependent ODE. Default is time dependent. +!~~~> Autonomous (1) or time dependent ODE (0). Default is time dependent. Autonomous = .NOT.(ICNTRL(1) == 0) !~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1) @@ -360,7 +360,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & AR_target_spc = ICNTRL(14) !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -528,7 +528,7 @@ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR) CASE (-6) PRINT * , '--> No of steps exceeds maximum bound' CASE (-7) - PRINT * , '--> Step size too small: T + 10*H = T', & + PRINT * , '--> Step size too small: T + 0.1*H = T', & ' or H < Roundoff' CASE (-8) PRINT * , '--> Matrix is repeatedly singular' @@ -587,9 +587,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -661,29 +658,28 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + & + ros_A((istage-1)*(istage-2)/2+j) * K(N*(j-1)+1:N*j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = Fcn(1:N) DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HC * K(N*(j-1)+1:N*j) END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HG * dFdT(1:N) END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) @@ -691,17 +687,15 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -717,7 +711,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -1018,13 +1011,11 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - call WCOPY(N,Fcn0,1,Fcn,1) - ! Fcn(1:N) = Fcn0(1:N) + Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ! K = 0.0_dp ! is this fix needed? hplin 14:04 -- not. 3 hours wiser later ELSEIF ( ros_NewF(istage) ) THEN - call WCOPY(N,Y,1,Ynew,1) - ! Ynew(1:N) = Y(1:N) + Ynew(1:N) = Y(1:N) DO j = 1, istage-1 ! In full vector space. Just use WAXPY as normal ! other entries in K are set to 1 previously. @@ -1090,7 +1081,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HC * K(N*(j-1)+SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HC * K(N*(j-1)+SPC_MAP(i)) ENDDO ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP) ! loop unrolling is consistently slower here. 18:58 @@ -1102,7 +1094,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HG * dFdT(SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HG * dFdT(SPC_MAP(i)) ENDDO ENDIF @@ -1153,8 +1146,7 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - !Y(1:N) = Ynew(1:N) + Y(1:N) = Ynew(1:N) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1482,13 +1474,11 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - call WCOPY(N,Fcn0,1,Fcn,1) - ! Fcn(1:N) = Fcn0(1:N) + Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ! K = 0.0_dp ! is this fix needed? hplin 14:04 -- not. 3 hours wiser later ELSEIF ( ros_NewF(istage) ) THEN - call WCOPY(N,Y,1,Ynew,1) - ! Ynew(1:N) = Y(1:N) + Ynew(1:N) = Y(1:N) DO j = 1, istage-1 ! In full vector space. Just use WAXPY as normal ! other entries in K are set to 1 previously. @@ -1554,7 +1544,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HC * K(N*(j-1)+SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HC * K(N*(j-1)+SPC_MAP(i)) ENDDO ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP) ! loop unrolling is consistently slower here. 18:58 @@ -1566,7 +1557,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) ! faster version: DO i = 1,rNVAR - K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) + HG * dFdT(SPC_MAP(i)) + K(ioffset+SPC_MAP(i)) = K(ioffset+SPC_MAP(i)) & + + HG * dFdT(SPC_MAP(i)) ENDDO ENDIF @@ -1616,8 +1608,7 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - !Y(1:N) = Ynew(1:N) + Y(1:N) = Ynew(1:N) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -1729,8 +1720,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1778,16 +1769,12 @@ SUBROUTINE ros_cPrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,rNVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:cNONZERO) = -Jac0(JVS_MAP(1:cNONZERO)) ghinv = ONE/(Direction*H*gam) DO i=1,rNVAR @@ -1873,8 +1860,6 @@ SUBROUTINE ros_cSolve( A, Pivot, b, map1, map2 ) Btmp = 0.d0 Atmp(map1(1:cNONZERO)) = A btmp(map2(1:rNVAR)) = b -! call cWCOPY(cNONZERO,LU_NONZERO,A,1,Atmp,1,map1) -! call cWCOPY(rNVAR,NVAR,B,1,Btmp,1,map2) CALL KppSolve( Atmp, btmp ) b = btmp(map2(1:rNVAR)) #endif @@ -1926,16 +1911,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_h211b_qssa.f90 b/int/rosenbrock_h211b_qssa.f90 index e925687..2d209ca 100644 --- a/int/rosenbrock_h211b_qssa.f90 +++ b/int/rosenbrock_h211b_qssa.f90 @@ -409,7 +409,7 @@ SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, & ENDIF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -662,9 +662,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp KPP_REAL, PARAMETER :: DeltaMin = 1.0E-5_dp -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -837,31 +834,25 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - !slim: CALL WCOPY(N,Fcn0,1,Fcn,1) Fcn(1:N) = Fcn0(1:N) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j = 1, istage-1 - ! CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - ! K(N*(j-1)+1),1,Ynew,1) - Ynew(1:N) = Ynew(1:N) + K(N*(j-1)+1:N*j) * ros_A((istage-1)*(istage-2)/2+j) + Ynew(1:N) = Ynew(1:N) + & + K(N*(j-1)+1:N*j) * ros_A((istage-1)*(istage-2)/2+j) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate( Tau, Ynew, Fcn ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 END IF ! if istage == 1 elseif ros_NewF(istage) - !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = Fcn(1:N) DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - ! CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + K(N*(j-1)+1:N*j) * HC END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) -! CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + dFdT(1:N) * HG END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) @@ -870,18 +861,14 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & !~~~> Compute the new solution - !slim: CALL WCOPY(N,Y,1,Ynew,1) Ynew(1:N) = Y(1:N) DO j=1,ros_S - ! CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) Ynew(1:N) = Ynew(1:N) + K(N*(j-1)+1:N*j) * ros_m(j) END DO !~~~> Compute the error estimation - !slim: CALL WSCAL(N,ZERO,Yerr,1) Yerr(1:N) = ZERO DO j=1,ros_S - ! CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) Yerr(1:N) = Yerr(1:N) + K(N*(j-1)+1:N*j) * ros_E(j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) @@ -908,7 +895,6 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, & ! new value is non-negative: Y = MAX(Ynew,ZERO) ELSE - !slim: CALL WCOPY(N,Ynew,1,Y,1) Y(1:N) = Ynew(1:N) ENDIF T = T + Direction*H @@ -995,9 +981,7 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - ! CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) dFdt(1:N) = dFdt(1:N) - ONE * FcN0(1:N) - ! CALL WSCAL(N,(ONE/Delta),dFdT,1) dFDT(1:N) = dFDT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -1045,16 +1029,12 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA - !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1) Ghimj = -Jac0 ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else - !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N diff --git a/int/rosenbrock_tlm.f90 b/int/rosenbrock_tlm.f90 index 4710de3..de5f3d4 100644 --- a/int/rosenbrock_tlm.f90 +++ b/int/rosenbrock_tlm.f90 @@ -372,7 +372,7 @@ SUBROUTINE RosenbrockTLM(N,Y,NTLM,Y_tlm, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -545,9 +545,6 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & LOGICAL :: RejectLastH, RejectMoreH, Singular !~~~> Local parameters KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 -!~~~> Locally called functions -! KPP_REAL WLAMCH -! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -625,21 +622,23 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & ioffset = N*(istage-1) ! Initialize stage solution - CALL WCOPY(N,Y,1,Ynew,1) - CALL WCOPY(N*NTLM,Y_tlm,1,Ynew_tlm,1) + Ynew(1:N) = Y(1:N) + Ynew_tlm(1:N,1:NTLM) = Y_tlm(1:N,1:NTLM) ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN - CALL WCOPY(N,Fcn0,1,Fcn,1) - CALL WCOPY(N*NTLM,Fcn0_tlm,1,Fcn_tlm,1) + Fcn(1:N) = Fcn0(1:N) + Fcn_tlm(1:N,1:NTLM) = Fcn0_tlm(1:N,1:NTLM) ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN DO j = 1, istage-1 - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K(N*(j-1)+1:N*j) DO itlm=1,NTLM - CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), & - K_tlm(N*(j-1)+1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Ynew_tlm(1:N,itlm) & + + ros_A((istage-1)*(istage-2)/2+j) & + * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO Tau = T + ros_Alpha(istage)*Direction*H @@ -651,29 +650,33 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & CALL Jac_SP_Vec ( Jac, Ynew_tlm(1,itlm), Fcn_tlm(1,itlm) ) END DO END IF ! if istage == 1 elseif ros_NewF(istage) - CALL WCOPY(N,Fcn,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = Fcn(1:N) DO itlm=1,NTLM - CALL WCOPY(N,Fcn_tlm(1,itlm),1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = Fcn_tlm(1:N,itlm) END DO DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) - CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) & + + HC * K(N*(j-1)+1:N*j) DO itlm=1,NTLM - CALL WAXPY(N,HC,K_tlm(N*(j-1)+1,itlm),1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + HC * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) - CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1) + K(ioffset+1:ioffset+N) = K(ioffset+1:ioffset+N) + HG * dFdT(1:N) DO itlm=1,NTLM CALL Jac_SP_Vec ( dJdT, Ynew_tlm(1,itlm), Tmp ) - CALL WAXPY(N,HG,Tmp,1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + HG * Tmp(1:N) END DO END IF CALL ros_Solve(Ghimj, Pivot, K(ioffset+1)) DO itlm=1,NTLM CALL Hess_Vec ( Hes0, K(ioffset+1), Y_tlm(1,itlm), Tmp ) - CALL WAXPY(N,ONE,Tmp,1,K_tlm(ioffset+1,itlm),1) + K_tlm(ioffset+1:ioffset+N,itlm) = K_tlm(ioffset+1:ioffset+N,itlm) & + + Tmp(1:N) CALL ros_Solve(Ghimj, Pivot, K_tlm(ioffset+1,itlm)) END DO @@ -681,29 +684,31 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & !~~~> Compute the new solution - CALL WCOPY(N,Y,1,Ynew,1) + Ynew(1:N) = Y(1:N) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1) + Ynew(1:N) = Ynew(1:N) + ros_M(j) * K(N*(j-1)+1:N*j) END DO DO itlm=1,NTLM - CALL WCOPY(N,Y_tlm(1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Y_tlm(1:N,itlm) DO j=1,ros_S - CALL WAXPY(N,ros_M(j),K_tlm(N*(j-1)+1,itlm),1,Ynew_tlm(1,itlm),1) + Ynew_tlm(1:N,itlm) = Ynew_tlm(1:N,itlm) & + + ros_M(j) * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO !~~~> Compute the error estimation - CALL Set2zero(N,Yerr) + Yerr(1:N) = 0.0_dp DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1) + Yerr(1:N) = Yerr(1:N) + ros_E(j) * K(N*(j-1)+1:N*j) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) IF (TLMtruncErr) THEN Err1 = 0.0d0 - CALL Set2zero(N*NTLM,Yerr_tlm) + Yerr_tlm(1:N,1:NTLM) = 0.0_dp DO itlm=1,NTLM DO j=1,ros_S - CALL WAXPY(N,ros_E(j),K_tlm(N*(j-1)+1,itlm),1,Yerr_tlm(1,itlm),1) + Yerr_tlm(1:N,itlm) = Yerr_tlm(1:N,itlm) & + + ros_E(j) * K_tlm(N*(j-1)+1:N*j,itlm) END DO END DO Err = ros_ErrorNorm_tlm(Y_tlm,Ynew_tlm,Yerr_tlm,AbsTol_tlm,RelTol_tlm,Err,VectorTol) @@ -717,8 +722,8 @@ SUBROUTINE ros_TLM_Int ( Y, NTLM, Y_tlm, & ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 - CALL WCOPY(N,Ynew,1,Y,1) - CALL WCOPY(N*NTLM,Ynew_tlm,1,Y_tlm,1) + Y(1:N) = Ynew(1:N) + Y_tlm(1:N,1:NTLM) = Ynew_tlm(1:N,1:NTLM) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step @@ -832,8 +837,8 @@ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate( T+Delta, Y, dFdT ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 - CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1) - CALL WSCAL(N,(ONE/Delta),dFdT,1) + dFdT(1:N) = dFdT(1:N) - Fcn0(1:N) + dFdT(1:N) = dFdT(1:N) * (ONE/Delta) END SUBROUTINE ros_FunTimeDerivative @@ -856,8 +861,8 @@ SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, Jac0, dJdT ) Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL JacTemplate( T+Delta, Y, dJdT ) ISTATUS(Njac) = ISTATUS(Njac) + 1 - CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1) - CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) - Jac0(1:LU_NONZERO) + dJdT(1:LU_NONZERO) = dJdT(1:LU_NONZERO) * (ONE/Delta) END SUBROUTINE ros_JacTimeDerivative @@ -895,8 +900,8 @@ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & DO WHILE (Singular) !~~~> Construct Ghimj = 1/(H*ham) - Jac0 - CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) - CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) + Ghimj(1:LU_NONZERO) = Jac0(1:LU_NONZERO) + Ghimj(1:LU_NONZERO) = -Ghimj(1:LU_NONZERO) ghinv = ONE/(Direction*H*gam) DO i=1,N Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv diff --git a/int/runge_kutta.f90 b/int/runge_kutta.f90 index 27bbab7..73f2eae 100644 --- a/int/runge_kutta.f90 +++ b/int/runge_kutta.f90 @@ -19,8 +19,9 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Precision USE KPP_ROOT_Parameters USE KPP_ROOT_Global - USE KPP_ROOT_Jacobian, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra + USE KPP_ROOT_Jacobian, ONLY : LU_DIAG + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppDecompCmplx, KppSolveCmplx IMPLICIT NONE PUBLIC @@ -391,7 +392,7 @@ SUBROUTINE RungeKutta( N,T,Tend,Y,RelTol,AbsTol, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -591,9 +592,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -640,9 +641,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -670,11 +671,11 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Prepare the loop-independent part of the right-hand side ! G = H*rkBgam(0)*F0 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - IF (rkMethod/=L3A) CALL WAXPY(N,rkBgam(0)*H, F0,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + IF (rkMethod/=L3A) G(1:N) = G(1:N) + rkBgam(0)*H * F0(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -683,12 +684,12 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -723,8 +724,8 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) - + Z4(1:N) = Z4(1:N) + DZ4(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT SDNewtonLoop @@ -742,21 +743,21 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN - CALL Set2Zero(N, DZ4) + DZ4(1:N) = 0.0_dp IF (rkMethod==L3A) THEN DZ4(1:N) = H*rkF(0)*F0(1:N) - IF (rkF(1) /= ZERO) CALL WAXPY(N, rkF(1), Z1, 1, DZ4, 1) - IF (rkF(2) /= ZERO) CALL WAXPY(N, rkF(2), Z2, 1, DZ4, 1) - IF (rkF(3) /= ZERO) CALL WAXPY(N, rkF(3), Z3, 1, DZ4, 1) + IF (rkF(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(1) * Z1(1:N) + IF (rkF(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(2) * Z2(1:N) + IF (rkF(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkF(3) * Z3(1:N) TMP = Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ1) - CALL WAXPY(N, H*rkBgam(4), DZ1, 1, DZ4, 1) + DZ4(1:N) = DZ4(1:N) + H*rkBgam(4) * DZ1(1:N) ELSE ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) END IF Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE @@ -790,9 +791,9 @@ SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum(d_i Z_i) - IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1) /= ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) @@ -978,33 +979,33 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,F0,Z1,Z2,Z3,R1,R2,R3) KPP_REAL :: T, H KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F0,F,R1,R2,R3,TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 IF (rkMethod==L3A) THEN - CALL WAXPY(N,-H*rkA(1,0),F0,1,R1,1) ! R1 <- R1 - h*A_10*F0 - CALL WAXPY(N,-H*rkA(2,0),F0,1,R2,1) ! R2 <- R2 - h*A_20*F0 - CALL WAXPY(N,-H*rkA(3,0),F0,1,R3,1) ! R3 <- R3 - h*A_30*F0 + R1(1:N) = R1(1:N) - H*rkA(1,0) * F0(1:N) ! R1 <- R1 - h*A_10*F0 + R2(1:N) = R2(1:N) - H*rkA(2,0) * F0(1:N) ! R2 <- R2 - h*A_20*F0 + R3(1:N) = R3(1:N) - H*rkA(3,0) * F0(1:N) ! R3 <- R3 - h*A_30*F0 END IF - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS diff --git a/int/runge_kutta_adj.f90 b/int/runge_kutta_adj.f90 index d7008ce..1074f8c 100644 --- a/int/runge_kutta_adj.f90 +++ b/int/runge_kutta_adj.f90 @@ -20,8 +20,10 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters USE KPP_ROOT_Global USE KPP_ROOT_Jacobian - USE KPP_ROOT_LinearAlgebra - + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppSolveTR, KppSolveCmplx, & + KppSolveTRCmplx, KppDecompCmplx, & + WGEFA, WGESL IMPLICIT NONE PUBLIC SAVE @@ -476,7 +478,7 @@ SUBROUTINE RungeKuttaADJ( N, Y, NADJ, Lambda,Tstart,Tend, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff = WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -810,15 +812,11 @@ SUBROUTINE rk_DPush( T, H, Y, Zstage, NewIt, E1, E2 )!, Jcb ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - ! CALL WCOPY(NVAR*3,Zstage,1,chk_Z(1,stack_ptr),1) chk_Y(1:N,stack_ptr) = Y(1:N) chk_Z(1:3*N,stack_ptr) = Zstage(1:3*N) chk_NiT( stack_ptr ) = NewIt IF (SaveLU) THEN #ifdef FULL_ALGEBRA - ! CALL WCOPY(NVAR*NVAR, E1,1,chk_E1(1,stack_ptr),1) - ! CALL WCOPYCmplx(NVAR*NVAR, E2,1,chk_E2(1,stack_ptr),1) DO j=1,NVAR DO i=1,NVAR chk_E1(NVAR*(j-1)+i,stack_ptr) = E1(i,j) @@ -826,13 +824,10 @@ SUBROUTINE rk_DPush( T, H, Y, Zstage, NewIt, E1, E2 )!, Jcb ) END DO END DO #else - ! CALL WCOPY(LU_NONZERO, E1,1,chk_E1(1,stack_ptr),1) - ! CALL WCOPYCmplx(LU_NONZERO, E2,1,chk_E2(1,stack_ptr),1) chk_E1(1:LU_NONZERO,stack_ptr) = E1(1:LU_NONZERO) chk_E2(1:LU_NONZERO,stack_ptr) = E2(1:LU_NONZERO) #endif END IF - !CALL WCOPY(LU_NONZERO,Jcb,1,chk_J(1,stack_ptr),1) END SUBROUTINE rk_DPush @@ -860,15 +855,11 @@ SUBROUTINE rk_DPop( T, H, Y, Zstage, NewIt, E1, E2 ) !, Jcb ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) - ! CALL WCOPY(NVAR*3,chk_Z(1,stack_ptr),1,Zstage,1) Zstage(1:3*NVAR) = chk_Z(1:3*NVAR,stack_ptr) NewIt = chk_NiT( stack_ptr ) IF (SaveLU) THEN #ifdef FULL_ALGEBRA - ! CALL WCOPY(NVAR*NVAR,chk_E1(1,stack_ptr),1, E1,1) - ! CALL WCOPYCmplx(NVAR*NVAR,chk_E2(1,stack_ptr),1, E2,1) DO j=1,NVAR DO i=1,NVAR E1(i,j) = chk_E1(NVAR*(j-1)+i,stack_ptr) @@ -876,13 +867,10 @@ SUBROUTINE rk_DPop( T, H, Y, Zstage, NewIt, E1, E2 ) !, Jcb ) END DO END DO #else - ! CALL WCOPY(LU_NONZERO,chk_E1(1,stack_ptr),1, E1,1) - ! CALL WCOPYCmplx(LU_NONZERO,chk_E2(1,stack_ptr),1, E2,1) E1(1:LU_NONZERO) = chk_E1(1:LU_NONZERO,stack_ptr) E2(1:LU_NONZERO) = chk_E2(1:LU_NONZERO,stack_ptr) #endif END IF - !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1) stack_ptr = stack_ptr - 1 @@ -903,9 +891,6 @@ SUBROUTINE rk_CPush(T, H, Y, dY, d2Y ) END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T - ! CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) - ! CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1) - ! CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1) chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR) chk_d2Y(1:NVAR,stack_ptr)= d2Y(1:NVAR) @@ -927,9 +912,6 @@ SUBROUTINE rk_CPop( T, H, Y, dY, d2Y ) END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) - ! CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) - ! CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1) - ! CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr) d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr) @@ -1034,9 +1016,9 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -1084,9 +1066,9 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -1122,11 +1104,11 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! G = H*rkBgam(0)*DZ4 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - CALL WAXPY(N,rkBgam(0)*H, DZ4,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + G(1:N) = G(1:N) + rkBgam(0)*H * DZ4(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -1135,12 +1117,13 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) + !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -1175,8 +1158,8 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) - + Z4(1:N) = Z4(1:N) + DZ4(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT SDNewtonLoop @@ -1196,12 +1179,12 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - CALL Set2Zero(N, DZ4) - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) - Err = RK_ErrorNorm(N,SCAL,DZ4) + DZ4(1:N) = 0.0_dp + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) + Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE CALL RK_ErrorEstimate(N,H,Y,T, & E1,IP1,Z1,Z2,Z3,SCAL,Err,FirstStep,Reject) @@ -1217,9 +1200,6 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ accept:IF (Err < ONE) THEN !~~~> STEP IS ACCEPTED IF (AdjointType == KPP_ROOT_discrete) THEN ! Save stage solution - ! CALL WCOPY(N,Z1,1,Zstage(1),1) - ! CALL WCOPY(N,Z2,1,Zstage(N+1),1) - ! CALL WCOPY(N,Z3,1,Zstage(2*N+1),1) Zstage(1:N) = Z1(1:N) Zstage(N+1:2*N) = Z2(1:N) Zstage(2*N+1:3*N) = Z3(1:N) @@ -1243,9 +1223,10 @@ SUBROUTINE RK_FwdIntegrator( N,Tstart,Tend,Y,AdjointType,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum (d_i Z_i) - IF (rkD(1) /= ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2) /= ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3) /= ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1) /= ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) + ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL) @@ -1367,12 +1348,12 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF !~~~> Jacobian values at stage vectors - CALL WADD(N,Y,Zstage(1),TMP) ! TMP <- Y + Z1 - CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) - CALL WADD(N,Y,Zstage(1+N),TMP) ! TMP <- Y + Z2 - CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) - CALL WADD(N,Y,Zstage(1+2*N),TMP) ! TMP <- Y + Z3 - CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) + TMP(1:N) = Y(1:N) + Zstage(1:N) ! TMP <- Y + Z1 + CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) + TMP(1:N) = Y(1:N) + Zstage(N+1:2*N) ! TMP <- Y + Z2 + CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) + TMP(1:N) = Y(1:N) + Zstage(2*N+1:3*N) ! TMP <- Y + Z3 + CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) END IF ! .not.Reject @@ -1398,7 +1379,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) Jbig(i,i) = Jbig(i,i) + ONE END DO CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,ISING) - ! CALL WGEFA(3*N,Jbig,IPbig,ISING) IF (ISING /= 0) THEN PRINT*,'Big guy is singular'; STOP END IF @@ -1435,7 +1415,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) DO i=1, 3*N Jbig(i,i) = ONE + Jbig(i,i) END DO - ! CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,ISING) CALL WGEFA(3*N,Jbig,IPbig,ISING) IF (ISING /= 0) THEN PRINT*,'Big guy is singular'; STOP @@ -1448,13 +1427,10 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Adj:DO iadj = 1, NADJ !~~~> Starting values for Newton iteration - ! CALL WCOPY(N,Lambda(1,iadj),1,U1(1,iadj),1) - ! CALL WCOPY(N,Lambda(1,iadj),1,U2(1,iadj),1) - ! CALL WCOPY(N,Lambda(1,iadj),1,U3(1,iadj),1) - CALL Set2Zero(N,U1(1,iadj)) - CALL Set2Zero(N,U2(1,iadj)) - CALL Set2Zero(N,U3(1,iadj)) - + U1(1:N,iadj) = 0.0_dp + U2(1:N,iadj) = 0.0_dp + U3(1:N,iadj) = 0.0_dp + !~~~> Initializations for Newton iteration NewtonDone = .FALSE. !~~~> Right Hand Side - part G for all Newton iterations @@ -1508,9 +1484,9 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF ! (AdjointSolve == Solve_adaptive) ! Update solution - CALL WAXPY(N,-ONE,DU1,1,U1(1,iadj),1) ! U1 <- U1 - DU1 - CALL WAXPY(N,-ONE,DU2,1,U2(1,iadj),1) ! U2 <- U2 - DU2 - CALL WAXPY(N,-ONE,DU3,1,U3(1,iadj),1) ! U3 <- U3 - DU3 + U1(1:N,iadj) = U1(1:N,iadj) - DU1(1:N) ! U1 <- U1 - DU1 + U2(1:N,iadj) = U2(1:N,iadj) - DU2(1:N) ! U2 <- U2 - DU2 + U3(1:N,iadj) = U3(1:N,iadj) - DU3(1:N) ! U3 <- U3 - DU3 IF (AdjointSolve == Solve_adaptive) THEN ! When performing an adaptive number of iterations @@ -1531,9 +1507,9 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) END IF ! Update adjoint solution: Y_adj <- Y_adj + sum (U_i) - CALL WAXPY(N,ONE,U1(1,iadj),1,Lambda(1,iadj),1) - CALL WAXPY(N,ONE,U2(1,iadj),1,Lambda(1,iadj),1) - CALL WAXPY(N,ONE,U3(1,iadj),1,Lambda(1,iadj),1) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U1(1:N,iadj) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U2(1:N,iadj) + Lambda(1:N,iadj) = Lambda(1:N,iadj) + U3(1:N,iadj) ELSE ! NewtonConverge = .false. @@ -1542,7 +1518,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) X(N+1:2*N) = -G2(1:N) X(2*N+1:3*N) = -G3(1:N) CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,ISING) - ! CALL WGESL('T',3*N,Jbig,IPbig,X) Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N) #else ! Commented lines for sparse big algebra: @@ -1555,7 +1530,6 @@ SUBROUTINE rk_DadjInt( NADJ,Lambda,Tstart,Tend,T,IERR ) X(1:N) = -G1(1:N) X(N+1:2*N) = -G2(1:N) X(2*N+1:3*N) = -G3(1:N) - ! CALL DGETRS('T',3*N,1,Jbig,3*N,IPbig,X,3*N,ISING) CALL WGESL('T',3*N,Jbig,IPbig,X) Lambda(1:N,iadj) = Lambda(1:N,iadj)+X(1:N)+X(N+1:2*N)+X(2*N+1:3*N) #endif @@ -1758,28 +1732,28 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,R1,R2,R3) KPP_REAL :: T,H KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F,R1,R2,R3,TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 - - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 - + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 + + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP, F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 + END SUBROUTINE RK_PrepareRHS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1802,46 +1776,48 @@ SUBROUTINE RK_PrepareRHS_Adj(N,H,Jac1,Jac2,Jac3,Lambda, & KPP_REAL, DIMENSION(N) :: F,TMP - CALL WCOPY(N,G1,1,R1,1) ! R1 <- G1 - CALL WCOPY(N,G2,1,R2,1) ! R2 <- G2 - CALL WCOPY(N,G3,1,R3,1) ! R3 <- G3 + R1(1:N) = G1(1:N) ! R1 <- G1 + R2(1:N) = G2(1:N) ! R2 <- G2 + R3(1:N) = G3(1:N) ! R3 <- G3 + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,1),U1,1,F,1) ! F1 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,1),U2,1,F,1) ! F1 <- F1 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,1),U3,1,F,1) ! F1 <- F1 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac1),F) #else - CALL JacTR_SP_Vec ( Jac1, F, TMP ) ! R1 <- -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + CALL JacTR_SP_Vec ( Jac1, F, TMP ) ! R1 <- -Jac(Y+Z1)^t*h*sum(A_j1*U_j) #endif - CALL WAXPY(N,ONE,U1,1,TMP,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) - CALL WAXPY(N,ONE,TMP,1,R1,1) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + TMP(1:N) = TMP(1:N) + U1(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R1(1:N) = R1(1:N) + TMP(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,2),U1,1,F,1) ! F2 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,2),U2,1,F,1) ! F2 <- F2 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,2),U3,1,F,1) ! F2 <- F2 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac2),F) #else - CALL JacTR_SP_Vec ( Jac2, F, TMP ) ! R2 <- -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + CALL JacTR_SP_Vec ( Jac2, F, TMP ) ! R2 <- -Jac(Y+Z2)^t*h*sum(A_j2*U_j) #endif - CALL WAXPY(N,ONE,U2,1,TMP,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) - CALL WAXPY(N,ONE,TMP,1,R2,1) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + TMP(1:N) = TMP(1:N) + U2(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R2(1:N) = R2(1:N) + TMP(1:N) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) + + F(1:N) = 0.0_dp + F(1:N) = F(1:N) - H*rkA(1,1) * U1(1:N) ! F1 <- -h*A_11*U1 + F(1:N) = F(1:N) - H*rkA(2,1) * U2(1:N) ! F1 <- F1 - h*A_21*U2 + F(1:N) = F(1:N) - H*rkA(3,1) * U3(1:N) ! F1 <- F1 - h*A_31*U3 - CALL SET2ZERO(N,F) - CALL WAXPY(N,-H*rkA(1,3),U1,1,F,1) ! F3 <- -h*A_11*U1 - CALL WAXPY(N,-H*rkA(2,3),U2,1,F,1) ! F3 <- F3 - h*A_21*U2 - CALL WAXPY(N,-H*rkA(3,3),U3,1,F,1) ! F3 <- F3 - h*A_31*U3 #ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac3),F) #else CALL JacTR_SP_Vec ( Jac3, F, TMP ) ! R3 <- -Jac(Y+Z3)^t*h*sum(A_j3*U_j) #endif - CALL WAXPY(N,ONE,U3,1,TMP,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j) - CALL WAXPY(N,ONE,TMP,1,R3,1) ! R3 <- U3 -Jac(Y+Z3)^t*h*sum(A_j3*U_j) - + TMP(1:N) = TMP(1:N) + U3(1:N) ! R1 <- U1 -Jac(Y+Z1)^t*h*sum(A_j1*U_j) + R3(1:N) = R3(1:N) + TMP(1:N) ! R2 <- U2 -Jac(Y+Z2)^t*h*sum(A_j2*U_j) END SUBROUTINE RK_PrepareRHS_Adj @@ -1864,29 +1840,29 @@ SUBROUTINE RK_PrepareRHS_G(N,H,Jac1,Jac2,Jac3,Lambda, & #endif KPP_REAL, DIMENSION(N) :: F - CALL SET2ZERO(N,G1) - CALL SET2ZERO(N,G2) - CALL SET2ZERO(N,G3) + G1(1:N) = 0.0_dp + G2(1:N) = 0.0_dp + G3(1:N) = 0.0_dp #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac1),Lambda) #else CALL JacTR_SP_Vec ( Jac1, Lambda, F ) ! F1 <- Jac(Y+Z1)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(1),F,1,G1,1) ! R1 <- R1 - h*B_1*F1 + G1(1:N) = G1(1:N) - H*rkB(1) * F(1:N) ! R1 <- R1 - h*B_1*F1 #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac2),Lambda) #else CALL JacTR_SP_Vec ( Jac2, Lambda, F ) ! F2 <- Jac(Y+Z2)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(2),F,1,G2,1) ! R2 <- R2 - h*B_2*F2 + G2(1:N) = G2(1:N) - H*rkB(2) * F(1:N) ! R2 <- R2 - h*B_2*F2 #ifdef FULL_ALGEBRA F = MATMUL(TRANSPOSE(Jac3),Lambda) #else CALL JacTR_SP_Vec ( Jac3, Lambda, F ) ! F3 <- Jac(Y+Z3)^t*Lambda #endif - CALL WAXPY(N,-H*rkB(3),F,1,G3,1) ! R3 <- R3 - h*B_3*F3 + G3(1:N) = G3(1:N) - H*rkB(3) * F(1:N) ! R3 <- R3 - h*B_3*F3 END SUBROUTINE RK_PrepareRHS_G diff --git a/int/runge_kutta_tlm.f90 b/int/runge_kutta_tlm.f90 index d0afc1f..b175444 100644 --- a/int/runge_kutta_tlm.f90 +++ b/int/runge_kutta_tlm.f90 @@ -21,6 +21,8 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Jacobian USE KPP_ROOT_LinearAlgebra + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, & + KppDecompCmplx, KppSolveCmplx IMPLICIT NONE PUBLIC @@ -436,7 +438,7 @@ SUBROUTINE RungeKuttaTLM( N, NTLM, Y, Y_tlm, T, Tend,RelTol,AbsTol, & END IF !~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0 - Roundoff=WLAMCH('E'); + Roundoff = EPSILON( 0.0_dp ) !~~~> Hmin = minimal step size IF (RCNTRL(1) == ZERO) THEN @@ -648,9 +650,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) !~~~> Starting values for Newton iteration IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN - CALL Set2zero(N,Z1) - CALL Set2zero(N,Z2) - CALL Set2zero(N,Z3) + Z1(1:N) = 0.0_dp + Z2(1:N) = 0.0_dp + Z3(1:N) = 0.0_dp ELSE ! Evaluate quadratic polynomial CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT) @@ -697,9 +699,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3 + Z1(1:N) = Z1(1:N) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2(1:N) = Z2(1:N) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3(1:N) = Z3(1:N) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -734,11 +736,11 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! G = H*rkBgam(0)*DZ4 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3 - CALL Set2Zero(N, G) - CALL WAXPY(N,rkBgam(0)*H, DZ4,1,G,1) - CALL WAXPY(N,rkTheta(1),Z1,1,G,1) - CALL WAXPY(N,rkTheta(2),Z2,1,G,1) - CALL WAXPY(N,rkTheta(3),Z3,1,G,1) + G(1:N) = 0.0_dp + G(1:N) = G(1:N) + rkBgam(0)*H * DZ4(1:N) + G(1:N) = G(1:N) + rkTheta(1) * Z1(1:N) + G(1:N) = G(1:N) + rkTheta(2) * Z2(1:N) + G(1:N) = G(1:N) + rkTheta(3) * Z3(1:N) !~~~> Initializations for Newton iteration NewtonDone = .FALSE. @@ -747,12 +749,13 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z4,TMP) ! TMP <- Y + Z4 + TMP(1:N) = Y(1:N) + Z4(1:N) ! TMP <- Y + Z4 CALL FUN_CHEM(T+H,TMP,DZ4) ! DZ4 <- Fun(Y+Z4) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N) - CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1) - CALL WAXPY (N, rkGamma/H, G,1, DZ4,1) + DZ4(1:N) = DZ4(1:N) - (rkGamma/H) * Z4(1:N) + DZ4(1:N) = DZ4(1:N) + (rkGamma/H) * G(1:N) + !~~~> Solve the linear system #ifdef FULL_ALGEBRA @@ -786,7 +789,7 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z4 <-- Z4 + DZ4 - CALL WAXPY(N,ONE,DZ4,1,Z4,1) + Z4(1:N) = Z4(1:N) + DZ4(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -807,11 +810,12 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (SdirkError) THEN ! DZ4(1:N) = rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4 - CALL Set2Zero(N, DZ4) - IF (rkD(1) /= ZERO) CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1) - IF (rkD(2) /= ZERO) CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1) - IF (rkD(3) /= ZERO) CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1) - CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1) + DZ4(1:N) = 0.0_dp + IF (rkD(1) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3) /= ZERO) DZ4(1:N) = DZ4(1:N) + rkD(3) * Z3(1:N) + DZ4(1:N) = DZ4(1:N) - Z4(1:N) + Err = RK_ErrorNorm(N,SCAL,DZ4) ELSE CALL RK_ErrorEstimate(N,H,Y,T, & @@ -824,11 +828,11 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) IF (Err < ONE) THEN !~~~> Jacobian values - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 CALL JAC_CHEM(T+rkC(1)*H,TMP,Jac1) ! Jac1 <- Jac(Y+Z1) - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 CALL JAC_CHEM(T+rkC(2)*H,TMP,Jac2) ! Jac2 <- Jac(Y+Z2) - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 CALL JAC_CHEM(T+rkC(3)*H,TMP,Jac3) ! Jac3 <- Jac(Y+Z3) TLMDIR: IF (TLMDirect) THEN @@ -853,7 +857,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Jbig(i,i) = ONE + Jbig(i,i) END DO !~~~> Solve the big system - ! CALL DGETRF(3*NVAR,3*NVAR,Jbig,3*NVAR,IPbig,j) CALL WGEFA(3*N,Jbig,IPbig,info) IF (info /= 0) THEN PRINT*,'Big big guy is singular'; STOP @@ -865,7 +868,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Zbig(2*NVAR+j) = Y_tlm(j,itlm) END DO Zbig = MATMUL(Ebig,Zbig) - !CALL DGETRS ('N',3*NVAR,1,Jbig,3*NVAR,IPbig,Zbig,3*NVAR,ISING) CALL WGESL('N',3*N,Jbig,IPbig,Zbig) DO j=1,NVAR Z1_tlm(j,itlm) = Zbig(j) @@ -907,7 +909,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) DO i=1, 3*N Jbig(i,i) = ONE + Jbig(i,i) END DO - ! CALL DGETRF(3*N,3*N,Jbig,3*N,IPbig,info) CALL WGEFA(3*N,Jbig,IPbig,info) IF (info /= 0) THEN PRINT*,'Big guy is singular'; STOP @@ -927,7 +928,6 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) ! Compute RHS CALL RK_PrepareRHS_TLMdirect(N,H,Jac1,Jac2,Jac3,Y_tlm(1,itlm),Zbig) ! Solve the system - ! CALL DGETRS('N',3*N,1,Jbig,3*N,IPbig,Zbig,3*N,ISING) CALL WGESL('N',3*N,Jbig,IPbig,Zbig) Z1_tlm(1:NVAR,itlm) = Zbig(1:NVAR) Z2_tlm(1:NVAR,itlm) = Zbig(NVAR+1:2*NVAR) @@ -944,9 +944,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Tlm:DO itlm = 1, NTLM !~~~> Starting values for Newton iteration - CALL Set2zero(N,Z1_tlm(1,itlm)) - CALL Set2zero(N,Z2_tlm(1,itlm)) - CALL Set2zero(N,Z3_tlm(1,itlm)) + Z1_tlm(1:N,itlm) = 0.0_dp + Z2_tlm(1:N,itlm) = 0.0_dp + Z3_tlm(1:N,itlm) = 0.0_dp !~~~> Initializations for Newton iteration IF (TLMNewtonEst) THEN @@ -998,9 +998,9 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) END IF !(TLMNewtonEst) ! Update solution - CALL WAXPY(N,-ONE,DZ1,1,Z1_tlm(1,itlm),1) ! Z1 <- Z1 - DZ1 - CALL WAXPY(N,-ONE,DZ2,1,Z2_tlm(1,itlm),1) ! Z2 <- Z2 - DZ2 - CALL WAXPY(N,-ONE,DZ3,1,Z3_tlm(1,itlm),1) ! Z3 <- Z3 - DZ3 + Z1_tlm(1:N,itlm) = Z1_tlm(1:N,itlm) - DZ1(1:N) ! Z1 <- Z1 - DZ1 + Z2_tlm(1:N,itlm) = Z2_tlm(1:N,itlm) - DZ2(1:N) ! Z2 <- Z2 - DZ2 + Z3_tlm(1:N,itlm) = Z3_tlm(1:N,itlm) - DZ3(1:N) ! Z3 <- Z3 - DZ3 ! Check error in Newton iterations IF (TLMNewtonEst) THEN @@ -1061,14 +1061,15 @@ SUBROUTINE RK_IntegratorTLM( N,NTLM,T,Tend,Y,Y_tlm,IERR ) Hold = H T = T+H ! Update solution: Y <- Y + sum(d_i Z_i) - IF (rkD(1)/=ZERO) CALL WAXPY(N,rkD(1),Z1,1,Y,1) - IF (rkD(2)/=ZERO) CALL WAXPY(N,rkD(2),Z2,1,Y,1) - IF (rkD(3)/=ZERO) CALL WAXPY(N,rkD(3),Z3,1,Y,1) + IF (rkD(1)/=ZERO) Y(1:N) = Y(1:N) + rkD(1) * Z1(1:N) + IF (rkD(2)/=ZERO) Y(1:N) = Y(1:N) + rkD(2) * Z2(1:N) + IF (rkD(3)/=ZERO) Y(1:N) = Y(1:N) + rkD(3) * Z3(1:N) + ! Update TLM solution: Y <- Y + sum(d_i*Z_i_tlm) DO itlm = 1,NTLM - IF (rkD(1)/=ZERO) CALL WAXPY(N,rkD(1),Z1_tlm(1,itlm),1,Y_tlm(1,itlm),1) - IF (rkD(2)/=ZERO) CALL WAXPY(N,rkD(2),Z2_tlm(1,itlm),1,Y_tlm(1,itlm),1) - IF (rkD(3)/=ZERO) CALL WAXPY(N,rkD(3),Z3_tlm(1,itlm),1,Y_tlm(1,itlm),1) + IF (rkD(1)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(1) * Z1_tlm(1:N,itlm) + IF (rkD(2)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(2) * Z2_tlm(1:N,itlm) + IF (rkD(3)/=ZERO) Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) + rkD(3) * Z3_tlm(1:N,itlm) END DO ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3 IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT) @@ -1264,27 +1265,27 @@ SUBROUTINE RK_PrepareRHS(N,T,H,Y,Z1,Z2,Z3,R1,R2,R3) KPP_REAL, INTENT(INOUT), DIMENSION(N) :: R1,R2,R3 KPP_REAL, DIMENSION(N) :: F, TMP - CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1 - CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2 - CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3 - - CALL WADD(N,Y,Z1,TMP) ! TMP <- Y + Z1 - CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 - - CALL WADD(N,Y,Z2,TMP) ! TMP <- Y + Z2 - CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 - - CALL WADD(N,Y,Z3,TMP) ! TMP <- Y + Z3 - CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + R1(1:N) = Z1(1:N) ! R1 <- Z1 + R2(1:N) = Z2(1:N) ! R2 <- Z2 + R3(1:N) = Z3(1:N) ! R3 <- Z3 + + TMP(1:N) = Y(1:N) + Z1(1:N) ! TMP <- Y + Z1 + CALL FUN_CHEM(T+rkC(1)*H,TMP,F) ! F1 <- Fun(Y+Z1) + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + + TMP(1:N) = Y(1:N) + Z2(1:N) ! TMP <- Y + Z2 + CALL FUN_CHEM(T+rkC(2)*H,TMP,F) ! F2 <- Fun(Y+Z2) + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 + + TMP(1:N) = Y(1:N) + Z3(1:N) ! TMP <- Y + Z3 + CALL FUN_CHEM(T+rkC(3)*H,TMP,F) ! F3 <- Fun(Y+Z3) + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS @@ -1308,39 +1309,40 @@ SUBROUTINE RK_PrepareRHS_TLM(N,H,Jac1,Jac2,Jac3,Y_tlm, & #endif KPP_REAL, DIMENSION(N) :: F, TMP - CALL WCOPY(N,Z1_tlm,1,R1,1) ! R1 <- Z1_tlm - CALL WCOPY(N,Z2_tlm,1,R2,1) ! R2 <- Z2_tlm - CALL WCOPY(N,Z3_tlm,1,R3,1) ! R3 <- Z3_tlm + R1(1:N) = Z1_tlm(1:N) ! R1 <- Z1_tlm + R2(1:N) = Z2_tlm(1:N) ! R2 <- Z2_tlm + R3(1:N) = Z3_tlm(1:N) ! R3 <- Z3_tlm - CALL WADD(N,Y_tlm,Z1_tlm,TMP) ! TMP <- Y + Z1 + TMP(1:N) = Y_tlm(1:N) + Z1_tlm(1:N) ! TMP <- Y + Z1 #ifdef FULL_ALGEBRA F = MATMUL(Jac1,TMP) #else - CALL Jac_SP_Vec ( Jac1, TMP, F ) ! F1 <- Jac(Y+Z1)*(Y_tlm+Z1_tlm) + CALL Jac_SP_Vec ( Jac1, TMP, F ) ! F1 <- Jac(Y+Z1)*(Y_tlm+Z1_tlm) #endif - CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1 - CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1 - CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1 + R1(1:N) = R1(1:N) - H*rkA(1,1) * F(1:N) ! R1 <- R1 - h*A_11*F1 + R2(1:N) = R2(1:N) - H*rkA(2,1) * F(1:N) ! R2 <- R2 - h*A_21*F1 + R3(1:N) = R3(1:N) - H*rkA(3,1) * F(1:N) ! R3 <- R3 - h*A_31*F1 + - CALL WADD(N,Y_tlm,Z2_tlm,TMP) ! TMP <- Y + Z2 + TMP(1:N) = Y_tlm(1:N) + Z2_tlm(1:N) ! TMP <- Y + Z2 #ifdef FULL_ALGEBRA F = MATMUL(Jac2,TMP) #else - CALL Jac_SP_Vec ( Jac2, TMP, F ) ! F2 <- Jac(Y+Z2)*(Y_tlm+Z2_tlm) + CALL Jac_SP_Vec ( Jac2, TMP, F ) ! F2 <- Jac(Y+Z2)*(Y_tlm+Z2_tlm) #endif - CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2 - CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2 - CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2 + R1(1:N) = R1(1:N) - H*rkA(1,2) * F(1:N) ! R1 <- R1 - h*A_12*F2 + R2(1:N) = R2(1:N) - H*rkA(2,2) * F(1:N) ! R2 <- R2 - h*A_22*F2 + R3(1:N) = R3(1:N) - H*rkA(3,2) * F(1:N) ! R3 <- R3 - h*A_32*F2 - CALL WADD(N,Y_tlm,Z3_tlm,TMP) ! TMP <- Y + Z3 + TMP(1:N) = Y_tlm(1:N) + Z3_tlm(1:N) ! TMP <- Y + Z3 #ifdef FULL_ALGEBRA F = MATMUL(Jac3,TMP) #else - CALL Jac_SP_Vec ( Jac3, TMP, F ) ! F3 <- Jac(Y+Z3)*(Y_tlm+Z3_tlm) + CALL Jac_SP_Vec ( Jac3, TMP, F ) ! F3 <- Jac(Y+Z3)*(Y_tlm+Z3_tlm) #endif - CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3 - CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3 - CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3 + R1(1:N) = R1(1:N) - H*rkA(1,3) * F(1:N) ! R1 <- R1 - h*A_13*F3 + R2(1:N) = R2(1:N) - H*rkA(2,3) * F(1:N) ! R2 <- R2 - h*A_23*F3 + R3(1:N) = R3(1:N) - H*rkA(3,3) * F(1:N) ! R3 <- R3 - h*A_33*F3 END SUBROUTINE RK_PrepareRHS_TLM diff --git a/int/sdirk.f90 b/int/sdirk.f90 index 490586f..9269541 100644 --- a/int/sdirk.f90 +++ b/int/sdirk.f90 @@ -22,9 +22,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WCOPY, WAXPY, & - WSCAL, WADD + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -366,7 +364,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -591,17 +589,16 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) - + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp !~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) IF (istage > 1) THEN DO j = 1, istage-1 ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) IF (StartNewton) THEN - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END IF END DO END IF @@ -613,13 +610,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) ! RHS <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N) - CALL WSCAL(N, H*rkGamma, RHS, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1) - CALL WAXPY (N, ONE, G,1, RHS,1) + RHS(1:N) = RHS(1:N) * (H * rkGamma) + RHS(1:N) = RHS(1:N) - Z(1:N,istage) + RHS(1:N) = RHS(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, RHS ) @@ -648,7 +645,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -675,9 +672,9 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF (sdMethod /= BEL) THEN ! All methods but Backward Euler - CALL Set2zero(N,TMP) + TMP(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) + IF (rkE(i)/=ZERO) TMP(1:N) = TMP(1:N) + rkE(i) * Z(1:N,i) END DO CALL SDIRK_Solve( H, N, E, IP, IER, TMP ) @@ -704,7 +701,7 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Ierr ) T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) DO i = 1,rkS - IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) + IF (rkD(i)/=ZERO) Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) END DO !~~~> Update scaling coefficients @@ -918,10 +915,10 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) KPP_REAL, INTENT(IN) :: E(LU_NONZERO) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) - KPP_REAL :: HGammaInv - HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + ! NOTE: This line reproduces the results of the + ! previous WAXPY call (@yantosca, 16 Oct 2025) + RHS(1:N) = RHS(1:N) * (ONE / (H * rkGamma)) #ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else diff --git a/int/sdirk4.f90 b/int/sdirk4.f90 index 7e1c04e..435cbe0 100644 --- a/int/sdirk4.f90 +++ b/int/sdirk4.f90 @@ -12,8 +12,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Global USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WAXPY, WCOPY + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -320,7 +319,7 @@ SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -553,11 +552,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & NewtonFactor(istage) = MAX(NewtonFactor(istage),Roundoff)**0.8d0 !~~~> STARTING VALUES FOR NEWTON ITERATION - CALL Set2zero(N,G) - CALL Set2zero(N,Z(1,istage)) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp IF (istage==1) THEN IF (FIRST.OR.NewtonReject) THEN - CALL Set2zero(N,Z(1,istage)) + Z(1:N,istage) = 0.0_dp ELSE W=ONE+rkGamma*H/Hold DO i=1,N @@ -567,12 +566,11 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & ELSE DO j = 1, istage-1 ! Gj(:) = sum_j Beta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkBeta(istage,j),Z(1,j),1,G,1) - ! CALL WAXPY(N,H*rkA(istage,j),FV(1,j),1,G,1) + G(1:N) = G(1:N) + rkBeta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END DO - IF (istage==5) CALL WCOPY(N,Z(1,istage),1,Yhat,1) ! Yhat(:) <- Z5(:) + IF (istage==5) Yhat(1:N) = Z(1:N,istage) ! Yhat(:) <- Z5(:) END IF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -597,7 +595,8 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & TMP(1:N) = Y(1:N) + Z(1:N,istage) CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) TMP(1:N) = G(1:N) - Z(1:N,istage) - CALL WAXPY(N,HGammaInv,TMP,1,RHS,1) ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:)) + ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:)) + RHS(1:N) = RHS(1:N) + HGammaInv * TMP(1:N) !~~~> SOLVE THE LINEAR SYSTEMS #ifdef FULL_ALGEBRA @@ -630,8 +629,9 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & END IF END IF NewtonErrOld = NewtonErr - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) ! Z(:) <-- Z(:)+RHS(:) - + ! Z(:) <-- Z(:)+RHS(:) + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) + END DO Newton !~~> END OF SIMPLIFIED NEWTON ITERATION @@ -671,13 +671,13 @@ SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, & Hold=H !~~~> COEFFICIENTS FOR CONTINUOUS SOLUTION - CALL WAXPY(N,ONE,Z(1,5),1,Y,1) ! Y(:) <-- Y(:)+Z5(:) - CALL WCOPY(N,Z(1,5),1,Yhat,1) ! Yhat <-- Z5 + Y(1:N) = Y(1:N) + Z(1:N,5) ! Y(:) <-- Y(:)+Z5(:) + Yhat(1:N) = Z(1:N,5) ! Yhat <-- Z5 DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj - CALL Set2zero(N,CONT(1,i)) + CONT(1:N,i) = 0.0_dp DO j = 1,5 - CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) + CONT(1:N,i) = CONT(1:N,i) + rkD(i,j) * Z(1:N,j) END DO END DO @@ -997,11 +997,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JS) - DO j=1,NVAR - DO j=1,NVAR - JV(i,j) = 0.0D0 - END DO - END DO + JV(1:NVAR,1:NVAR) = 0.0d0 DO i=1,LU_NONZERO JV(LU_IROW(i),LU_ICOL(i)) = JS(i) END DO diff --git a/int/sdirk_adj.f90 b/int/sdirk_adj.f90 index 8f58fa0..59f9b53 100644 --- a/int/sdirk_adj.f90 +++ b/int/sdirk_adj.f90 @@ -20,14 +20,13 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG USE KPP_ROOT_Jacobian, ONLY: Jac_SP_Vec, JacTR_SP_Vec - USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, & - KppSolveTR, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD - + USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, KppSolveTR + IMPLICIT NONE PUBLIC SAVE -!~~~> Flags to determine if we should call the UPDATE_* routines from within +!~~~> Flags to determine if we should call the UPDATE_* routines from within !~~~> the integrator. If using KPP in an external model, you might want to !~~~> disable these calls (via ICNTRL(15)) to avoid excess computations. LOGICAL, PRIVATE :: Do_Update_RCONST @@ -38,7 +37,7 @@ MODULE KPP_ROOT_Integrator INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, & Nrej=5, Ndec=6, Nsol=7, Nsng=8, & Ntexit=1, Nhexit=2, Nhnew=3 - + CONTAINS SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & @@ -56,7 +55,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & !~~~> NADJ - No. of cost functionals for which adjoints ! are evaluated simultaneously ! If single cost functional is considered (like in -! most applications) simply set NADJ = 1 +! most applications) simply set NADJ = 1 INTEGER :: NADJ !~~~> Lambda - Sensitivities w.r.t. concentrations ! Note: Lambda (1:NVAR,j) contains sensitivities of @@ -64,7 +63,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & KPP_REAL :: Lambda(NVAR,NADJ) !~~~> Tolerances for adjoint calculations ! (used for full continuous adjoint, and for controlling -! iterations when used to solve the discrete adjoint) +! iterations when used to solve the discrete adjoint) KPP_REAL, INTENT(IN) :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ) KPP_REAL, INTENT(IN) :: TIN ! Start Time KPP_REAL, INTENT(IN) :: TOUT ! End Time @@ -89,7 +88,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & ICNTRL(5) = 8 ! Max no. of Newton iterations ICNTRL(7) = 1 ! Adjoint solution by: 0=Newton, 1=direct ICNTRL(8) = 1 ! Save fwd LU factorization: 0 = do *not* save, 1 = save - ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. + ICNTRL(15) = 5 ! Call Update_SUN and Update_RCONST from w/in the int. !~~~> if optional parameters are given, and if they are /= 0, !~~~> then use them to overwrite default settings @@ -99,7 +98,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & IF ( PRESENT( RCNTRL_U ) ) THEN WHERE( RCNTRL_U > 0 ) RCNTRL = RCNTRL_U ENDIF - + !~~~> Determine the settings of the Do_Update_* flags, which determine !~~~> whether or not we need to call Update_* routines in the integrator !~~~> (or not, if we are calling them from a higher-level) @@ -109,7 +108,7 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & ! = 2 ! Call Update_PHOTO from within the integrator ! = 3 ! Call Update_RCONST and Update_PHOTO from w/in the int. ! = 4 ! Call Update_SUN from within the integrator - ! = 5 ! Call Update_SUN and Update_RCONST from within the int. + ! = 5 ! Call Update_SUN and Update_RCONST from within the int. ! = 6 ! Call Update_SUN and Update_PHOTO from within the int. ! = 7 ! Call Update_SUN, Update_PHOTO, Update_RCONST w/in int. CALL Integrator_Update_Options( ICNTRL(15), & @@ -136,12 +135,12 @@ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & !~~~> Debug option: number of steps ! Ntotal = Ntotal + ISTATUS(Nstp) ! WRITE(6,777) ISTATUS(Nstp),Ntotal,VAR(ind_O3),VAR(ind_NO2) - ! 777 FORMAT('NSTEPS=',I5,' (',I5,') O3=',E24.14,' NO2=',E24.14) + ! 777 FORMAT('NSTEPS=',I5,' (',I5,') O3=',E24.14,' NO2=',E24.14) IF (Ierr < 0) THEN PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (Ierr=',Ierr,')' ENDIF - + ! if optional parameters are given for output ! use them to store information in them IF ( PRESENT( ISTATUS_U ) ) ISTATUS_U = ISTATUS @@ -168,14 +167,14 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ! This code is based on the SDIRK4 routine in the above book. ! ! Methods: -! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 -! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant -! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 +! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 +! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant +! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! ! (C) Adrian Sandu, July 2005 ! Virginia Polytechnic Institute and State University ! Contact: sandu@cs.vt.edu -! Revised by Philipp Miehe and Adrian Sandu, May 2006 +! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! This implementation is part of KPP - the Kinetic PreProcessor !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! @@ -210,7 +209,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. -!~~~> +!~~~> ! ICNTRL(1) = not used ! ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors @@ -299,16 +298,16 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE -! Arguments +! Arguments INTEGER, INTENT(IN) :: N, NADJ, ICNTRL(20) KPP_REAL, INTENT(INOUT) :: Y(NVAR), Lambda(NVAR,NADJ) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, & RelTol(NVAR), AbsTol(NVAR), RCNTRL(20), & RelTol_adj(NVAR,NADJ), AbsTol_adj(NVAR,NADJ) INTEGER, INTENT(OUT) :: Ierr - INTEGER, INTENT(INOUT) :: ISTATUS(20) + INTEGER, INTENT(INOUT) :: ISTATUS(20) KPP_REAL, INTENT(OUT) :: RSTATUS(20) - + !~~~> SDIRK method coefficients, up to 5 stages INTEGER, PARAMETER :: Smax = 5 INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5 @@ -328,16 +327,16 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & #else KPP_REAL, DIMENSION(:,:), POINTER :: chk_J #endif -! Local variables +! Local variables KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, & FacMin, Facmax, FacSafe, FacRej, & ThetaMin, NewtonTol, Qmin, Qmax - LOGICAL :: SaveLU, DirectADJ + LOGICAL :: SaveLU, DirectADJ INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 - + !~~~> Initialize statistics ISTATUS(1:20) = 0 RSTATUS(1:20) = ZERO @@ -351,7 +350,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & ITOL = 0 END IF -!~~~> ICNTRL(3) - method selection +!~~~> ICNTRL(3) - method selection SELECT CASE (ICNTRL(3)) CASE (0,1) CALL Sdirk2a @@ -366,7 +365,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & CASE DEFAULT CALL Sdirk2a END SELECT - + !~~~> The maximum number of time steps admitted IF (ICNTRL(4) == 0) THEN Max_no_steps = 200000 @@ -376,7 +375,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4) CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr) END IF - + !~~~> The maximum number of Newton iterations admitted IF(ICNTRL(5) == 0)THEN NewtonMaxit=8 @@ -388,14 +387,14 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & END IF END IF -!~~~> Solve ADJ equations directly or by Newton iterations +!~~~> Solve ADJ equations directly or by Newton iterations DirectADJ = (ICNTRL(7) == 1) - + !~~~> Save or not the forward LU factorization SaveLU = (ICNTRL(8) /= 0) .AND. (.NOT.DirectADJ) !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -406,7 +405,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Upper bound on the step size: (positive value) IF (RCNTRL(2) == ZERO) THEN Hmax = ABS(Tfinal-Tinitial) @@ -416,7 +415,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN Hstart = MAX(Hmin,Roundoff) @@ -426,7 +425,7 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF - + !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax IF (RCNTRL(4) == ZERO) THEN FacMin = 0.2_dp @@ -505,19 +504,19 @@ SUBROUTINE SDIRKADJ(N, NADJ, Tinitial, Tfinal, Y, Lambda, & END IF END DO END IF - + IF (Ierr < 0) RETURN - -!~~~> Allocate memory buffers + +!~~~> Allocate memory buffers CALL SDIRK_AllocBuffers -!~~~> Call forward integration +!~~~> Call forward integration CALL SDIRK_FwdInt( N, Tinitial, Tfinal, Y, Ierr ) -!~~~> Call adjoint integration +!~~~> Call adjoint integration CALL SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) -!~~~> Free memory buffers +!~~~> Free memory buffers CALL SDIRK_FreeBuffers @@ -533,26 +532,26 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) USE KPP_ROOT_Parameters IMPLICIT NONE -!~~~> Arguments: +!~~~> Arguments: INTEGER :: N KPP_REAL, INTENT(INOUT) :: Y(NVAR) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: - KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & + +!~~~> Local variables: + KPP_REAL :: Z(NVAR,Smax), G(NVAR), TMP(NVAR), & NewtonRate, SCAL(NVAR), RHS(NVAR), & T, H, Theta, Hratio, NewtonPredictedErr, & Qnewton, Err, Fac, Hnew,Tdirection, & NewtonIncrement, NewtonIncrementOld INTEGER :: j, ISING, istage, NewtonIter, IP(NVAR) LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA + +#ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: FJAC, E -#else +#else KPP_REAL, DIMENSION(LU_NONZERO):: FJAC, E -#endif +#endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -585,14 +584,14 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) IF (ISING /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF - END IF + END IF IF (ISTATUS(Nstp) > Max_no_steps) THEN CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN - END IF + END IF IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN - END IF + END IF stages:DO istage = 1, rkS @@ -601,47 +600,48 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) - -!~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp + +!~~~> Prepare the loop-independent part of the right-hand side IF (istage > 1) THEN - DO j = 1, istage-1 - ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) - ! Zi(:) = sum_j Alpha(i,j)*Zj(:) - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) - END DO + DO j = 1, istage-1 + ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) + + ! Zi(:) = sum_j Alpha(i,j)*Zj(:) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) + END DO END IF !~~~> Initializations for Newton iteration NewtonDone = .FALSE. Fac = 0.5d0 ! Step reduction factor if too many iterations - + NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS) ! RHS <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! RHS(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*RHS(1:N) - CALL WSCAL(N, H*rkGamma, RHS, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, RHS, 1) - CALL WAXPY (N, ONE, G,1, RHS,1) + RHS(1:N) = RHS(1:N) * (H * rkGamma) + RHS(1:N) = RHS(1:N) - Z(1:N,istage) + RHS(1:N) = RHS(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( 'N', H, N, E, IP, ISING, RHS ) - + !~~~> Check convergence of Newton iterations CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 + NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process + ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) IF (NewtonPredictedErr >= NewtonTol) THEN @@ -656,14 +656,14 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+RHS(:) - CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) - + Z(1:N,istage) = Z(1:N,istage) + RHS(1:N) + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT NewtonLoop - + END DO NewtonLoop - + IF (.NOT.NewtonDone) THEN !CALL RK_ErrorMsg(-12,T,H,Ierr); H = Fac*H; Reject=.TRUE. @@ -673,7 +673,7 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> End of implified Newton iterations - + END DO stages @@ -681,10 +681,10 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - CALL Set2zero(N,TMP) + TMP(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,TMP,1) - END DO + IF (rkE(i)/=ZERO) TMP(1:N) = TMP(1:N) + rkE(i) * Z(1:N,i) + END DO CALL SDIRK_Solve( 'N', H, N, E, IP, ISING, TMP ) CALL SDIRK_ErrorNorm(N, TMP, SCAL, Err) @@ -708,10 +708,10 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) !~~~> Update time and solution T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO i = 1,rkS - IF (rkD(i)/=ZERO) CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) - END DO - + DO i = 1,rkS + IF (rkD(i)/=ZERO) Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) + END DO + !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL) @@ -746,16 +746,16 @@ SUBROUTINE SDIRK_FwdInt( N,Tinitial,Tfinal,Y,Ierr ) END IF Reject = .TRUE. SkipJac = .TRUE. - SkipLU = .FALSE. + SkipLU = .FALSE. IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1 - + END IF accept - + END DO Tloop ! Successful return Ierr = 1 - + END SUBROUTINE SDIRK_FwdInt @@ -767,34 +767,34 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) USE KPP_ROOT_Parameters IMPLICIT NONE -!~~~> Arguments: +!~~~> Arguments: INTEGER, INTENT(IN) :: N, NADJ KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ) INTEGER, INTENT(OUT) :: Ierr - -!~~~> Local variables: + +!~~~> Local variables: KPP_REAL :: Y(NVAR) KPP_REAL :: Z(NVAR,Smax), U(NVAR,NADJ,Smax), & - TMP(NVAR), G(NVAR), & + TMP(NVAR), G(NVAR), & NewtonRate, SCAL(NVAR), DU(NVAR), & T, H, Theta, NewtonPredictedErr, & NewtonIncrement, NewtonIncrementOld INTEGER :: j, ISING, istage, iadj, NewtonIter, & IP(NVAR), IP_adj(NVAR) LOGICAL :: Reject, SkipJac, SkipLU, NewtonDone - -#ifdef FULL_ALGEBRA + +#ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: E, Jac, E_adj -#else +#else KPP_REAL, DIMENSION(LU_NONZERO):: E, Jac, E_adj -#endif +#endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Time loop begins !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Tloop: DO WHILE ( stack_ptr > 0 ) - + !~~~> Recover checkpoints for stage values and vectors CALL SDIRK_Pop( T, H, Y, Z, E, IP ) @@ -806,7 +806,7 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) IF (ISING /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF - END IF + END IF stages:DO istage = rkS, 1, -1 @@ -814,8 +814,8 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) TMP(1:N) = Y(1:N) + Z(1:N,istage) CALL JAC_CHEM(T+rkC(istage)*H,TMP,Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 - - IF (DirectADJ) THEN + + IF (DirectADJ) THEN #ifdef FULL_ALGEBRA E_adj(1:N,1:N) = -Jac(1:N,1:N) DO j=1,N @@ -838,61 +838,61 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) END IF adj: DO iadj = 1, NADJ - + !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(ITOL, AbsTol_adj(1:NVAR,iadj), & RelTol_adj(1:NVAR,iadj), Lambda(1:NVAR,iadj), SCAL) - + !~~~> Prepare the loop-independent part of the right-hand side ! G(:) = H*Jac^T*( B(i)*Lambda + sum_j A(j,i)*Uj(:) ) G(1:N) = rkB(istage)*Lambda(1:N,iadj) IF (istage < rkS) THEN DO j = istage+1, rkS - CALL WAXPY(N,rkA(j,istage),U(1,iadj,j),1,G,1) + G(1:N) = G(1:N) + rkA(j,istage) * U(1:N,iadj,j) END DO END IF -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA TMP = MATMUL(TRANSPOSE(Jac),G) ! DZ <- Jac(Y+Z)*Y_tlm -#else - CALL JacTR_SP_Vec ( Jac, G, TMP ) -#endif +#else + CALL JacTR_SP_Vec ( Jac, G, TMP ) +#endif G(1:N) = H*TMP(1:N) -DirADJ:IF (DirectADJ) THEN +DirADJ:IF (DirectADJ) THEN CALL SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, ISING, G ) U(1:N,iadj,istage) = G(1:N) - + ELSE DirADJ !~~~> Initializations for Newton iteration - CALL Set2zero(N,U(1,iadj,istage)) + U(1:N,iadj,istage) = 0.0_dp NewtonDone = .FALSE. - + NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side -#ifdef FULL_ALGEBRA - TMP = MATMUL(TRANSPOSE(Jac),U(1:N,iadj,istage)) -#else - CALL JacTR_SP_Vec ( Jac, U(1:N,iadj,istage), TMP ) -#endif +#ifdef FULL_ALGEBRA + TMP = MATMUL(TRANSPOSE(Jac),U(1:N,iadj,istage)) +#else + CALL JacTR_SP_Vec ( Jac, U(1:N,iadj,istage), TMP ) +#endif DU(1:N) = U(1:N,iadj,istage) - (H*rkGamma)*TMP(1:N) - G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( 'T', H, N, E, IP, ISING, DU ) - + !~~~> Check convergence of Newton iterations - + CALL SDIRK_ErrorNorm(N, DU, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) - NewtonRate = 2.0d0 + NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) - ! Predict error at the end of Newton process + ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) ! Non-convergence of Newton: predicted error too large @@ -904,18 +904,18 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) NewtonIncrementOld = NewtonIncrement ! Update solution U(1:N,iadj,istage) = U(1:N,iadj,istage) - DU(1:N) - + ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) ! AbsTol is often inappropriate for adjoints - ! we do at least 4 Newton iterations to ensure convergence ! of all adjoint components IF ((NewtonIter>=4) .AND. NewtonDone) EXIT NewtonLoop - + END DO NewtonLoop - + !~~~> If Newton iterations fail employ the direct solution - IF (.NOT.NewtonDone) THEN + IF (.NOT.NewtonDone) THEN PRINT*,'Problems with Newton Adjoint!!!' #ifdef FULL_ALGEBRA E_adj(1:N,1:N) = -Jac(1:N,1:N) @@ -938,31 +938,30 @@ SUBROUTINE SDIRK_DadjInt( N, NADJ, Lambda, Ierr ) END IF CALL SDIRK_Solve ( 'T', H, N, E_adj, IP_adj, ISING, G ) U(1:N,iadj,istage) = G(1:N) - + END IF !~~~> End of simplified Newton iterations - + END IF DirADJ END DO adj - + END DO stages !~~~> Update adjoint solution ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) - DO istage = 1,rkS + DO istage = 1,rkS DO iadj = 1,NADJ Lambda(1:N,iadj) = Lambda(1:N,iadj) + U(1:N,iadj,istage) - !CALL WAXPY(N,ONE,U(1:N,iadj,istage),1,Lambda(1,iadj),1) - END DO - END DO + END DO + END DO END DO Tloop ! Successful return Ierr = 1 - + END SUBROUTINE SDIRK_DadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -970,23 +969,23 @@ SUBROUTINE SDIRK_AllocBuffers !~~~> Allocate buffer space for checkpointing !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i - + ALLOCATE( chk_H(Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer H'; STOP - END IF + END IF ALLOCATE( chk_T(Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer T'; STOP - END IF + END IF ALLOCATE( chk_Y(NVAR,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer Y'; STOP - END IF + END IF ALLOCATE( chk_Z(NVAR,rkS,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer K'; STOP - END IF + END IF IF (SaveLU) THEN #ifdef FULL_ALGEBRA ALLOCATE( chk_J(NVAR,NVAR,Max_no_steps), STAT=i ) @@ -995,13 +994,13 @@ SUBROUTINE SDIRK_AllocBuffers #endif IF (i/=0) THEN PRINT*,'Failed allocation of buffer J'; STOP - END IF + END IF ALLOCATE( chk_P(NVAR,Max_no_steps), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer P'; STOP - END IF - END IF - + END IF + END IF + END SUBROUTINE SDIRK_AllocBuffers @@ -1010,34 +1009,34 @@ SUBROUTINE SDIRK_FreeBuffers !~~~> Dallocate buffer space for discrete adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i - + DEALLOCATE( chk_H, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer H'; STOP - END IF + END IF DEALLOCATE( chk_T, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer T'; STOP - END IF + END IF DEALLOCATE( chk_Y, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer Y'; STOP - END IF + END IF DEALLOCATE( chk_Z, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer K'; STOP - END IF + END IF IF (SaveLU) THEN DEALLOCATE( chk_J, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer J'; STOP - END IF + END IF DEALLOCATE( chk_P, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer P'; STOP - END IF - END IF - + END IF + END IF + END SUBROUTINE SDIRK_FreeBuffers @@ -1047,35 +1046,35 @@ SUBROUTINE SDIRK_Push( T, H, Y, Z, E, P ) !~~~> Saves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - KPP_REAL :: T, H, Y(NVAR), Z(NVAR,Smax) + KPP_REAL :: T, H, Y(NVAR), Z(NVAR,Smax) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) -#else +#else KPP_REAL :: E(LU_NONZERO) #endif - + stack_ptr = stack_ptr + 1 IF ( stack_ptr > Max_no_steps ) THEN PRINT*,'Push failed: buffer overflow' STOP - END IF + END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_Z(1:NVAR,1:rkS,stack_ptr) = Z(1:NVAR,1:rkS) - IF (SaveLU) THEN + IF (SaveLU) THEN #ifdef FULL_ALGEBRA chk_J(1:NVAR,1:NVAR,stack_ptr) = E(1:NVAR,1:NVAR) chk_P(1:NVAR,stack_ptr) = P(1:NVAR) -#else +#else chk_J(1:LU_NONZERO,stack_ptr) = E(1:LU_NONZERO) #endif END IF - + END SUBROUTINE SDIRK_Push - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1086,14 +1085,14 @@ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) -#else +#else KPP_REAL :: E(LU_NONZERO) #endif - + IF ( stack_ptr <= 0 ) THEN PRINT*,'Pop failed: empty buffer' STOP - END IF + END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) @@ -1102,13 +1101,13 @@ SUBROUTINE SDIRK_Pop( T, H, Y, Z, E, P ) #ifdef FULL_ALGEBRA E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr) P(1:NVAR) = chk_P(1:NVAR,stack_ptr) -#else +#else E(1:LU_NONZERO) = chk_J(1:LU_NONZERO,stack_ptr) #endif END IF stack_ptr = stack_ptr - 1 - + END SUBROUTINE SDIRK_Pop !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1128,14 +1127,14 @@ SUBROUTINE SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL) END DO END IF END SUBROUTINE SDIRK_ErrorScale - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! +! INTEGER :: i, N - KPP_REAL :: Y(N), SCAL(N), Err + KPP_REAL :: Y(N), SCAL(N), Err Err = ZERO DO i=1,N Err = Err+(Y(i)*SCAL(i))**2 @@ -1145,7 +1144,7 @@ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) END SUBROUTINE SDIRK_ErrorNorm -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Handles all error messages @@ -1183,17 +1182,17 @@ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) PRINT *, "T=", T, "and H=", H END SUBROUTINE SDIRK_ErrorMsg - - + + !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipJac, SkipLU, E, IP, Reject, ISING ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + IMPLICIT NONE - + KPP_REAL, INTENT(INOUT) :: H KPP_REAL, INTENT(IN) :: T, Y(NVAR) LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject @@ -1210,9 +1209,9 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & ConsecutiveSng = 0 ISING = 1 - + Hloop: DO WHILE (ISING /= 0) - + HGammaInv = ONE/(H*rkGamma) !~~~> Compute the Jacobian @@ -1222,8 +1221,8 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & IF (.NOT. SkipJac) THEN CALL JAC_CHEM( T, Y, FJAC ) ISTATUS(Njac) = ISTATUS(Njac) + 1 - END IF - + END IF + #ifdef FULL_ALGEBRA DO j=1,NVAR DO i=1,NVAR @@ -1253,7 +1252,7 @@ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipLU = .FALSE. Reject = .TRUE. END IF - + END DO Hloop END SUBROUTINE SDIRK_PrepareMatrix @@ -1276,18 +1275,18 @@ SUBROUTINE SDIRK_Solve ( Transp, H, N, E, IP, ISING, RHS ) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) KPP_REAL :: HGammaInv - + HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + RHS(1:N) = RHS(1:N) * HGammaInv SELECT CASE (TRANSP) CASE ('N') -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolve(E, RHS) #endif CASE ('T') -#ifdef FULL_ALGEBRA +#ifdef FULL_ALGEBRA CALL DGETRS( 'T', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolveTR(E, RHS, RHS) @@ -1297,7 +1296,7 @@ SUBROUTINE SDIRK_Solve ( Transp, H, N, E, IP, ISING, RHS ) STOP END SELECT ISTATUS(Nsol) = ISTATUS(Nsol) + 1 - + END SUBROUTINE SDIRK_Solve @@ -1385,7 +1384,7 @@ SUBROUTINE Sdirk4a rkAlpha(5,2) = 6.559571569643355712998131800797873d0 rkAlpha(5,3) = -15.90772144271326504260996815012482d0 rkAlpha(5,4) = 25.34908987169226073668861694892683d0 - + !~~~> Coefficients for continuous solution ! rkD(1,1)= 24.74416644927758d0 ! rkD(1,2)= -4.325375951824688d0 @@ -1414,9 +1413,9 @@ SUBROUTINE Sdirk4a ! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) ! END DO ! END DO - + rkELO = 4.0d0 - + END SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1503,7 +1502,7 @@ SUBROUTINE Sdirk4b rkAlpha(5,2) = 9.d0 rkAlpha(5,3) = -56.81818181818181818181818181818182d0 rkAlpha(5,4) = 54.d0 - + END SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1545,7 +1544,7 @@ SUBROUTINE Sdirk2a ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 3.414213562373095048801688724209698d0 - + END SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1587,7 +1586,7 @@ SUBROUTINE Sdirk2b ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = .5857864376269049511983112757903019d0 - + END SUBROUTINE Sdirk2b @@ -1661,16 +1660,16 @@ SUBROUTINE FUN_CHEM( T, Y, P ) KPP_REAL :: T, Told KPP_REAL :: Y(NVAR), P(NVAR) - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() IF ( Do_Update_RCONST ) CALL Update_RCONST() - + CALL Fun( Y, FIX, RCONST, P ) - + TIME = Told - + END SUBROUTINE FUN_CHEM @@ -1691,7 +1690,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #else KPP_REAL :: JV(LU_NONZERO) #endif - + Told = TIME TIME = T IF ( Do_Update_SUN ) CALL Update_SUN() @@ -1699,11 +1698,7 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JS) - DO j=1,NVAR - DO j=1,NVAR - JV(i,j) = 0.0d0 - END DO - END DO + JV(1:NVAR,1:NVAR) = 0.0d0 DO i=1,LU_NONZERO JV(LU_IROW(i),LU_ICOL(i)) = JS(i) END DO @@ -1715,5 +1710,3 @@ SUBROUTINE JAC_CHEM( T, Y, JV ) END SUBROUTINE JAC_CHEM END MODULE KPP_ROOT_Integrator - - diff --git a/int/sdirk_tlm.f90 b/int/sdirk_tlm.f90 index 92d161f..1605570 100644 --- a/int/sdirk_tlm.f90 +++ b/int/sdirk_tlm.f90 @@ -20,9 +20,7 @@ MODULE KPP_ROOT_Integrator USE KPP_ROOT_Parameters USE KPP_ROOT_JacobianSP, ONLY : LU_DIAG USE KPP_ROOT_Jacobian, ONLY : Jac_SP_Vec - USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve, Set2zero, & - WLAMCH, WCOPY, WAXPY, & - WSCAL, WADD + USE KPP_ROOT_LinearAlgebra, ONLY : KppDecomp, KppSolve IMPLICIT NONE PUBLIC @@ -396,7 +394,7 @@ SUBROUTINE SdirkTLM(N, NTLM, Tinitial, Tfinal, Y, Y_tlm, RelTol, AbsTol, & END IF !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN @@ -596,17 +594,17 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z(1,istage)) + G(1:N) = 0.0_dp + Z(1:N,istage) = 0.0_dp -!~~~> Prepare the loop-independent part of the right-hand side - CALL Set2zero(N,G) +!~~~> Prepare the loop-independent part of the right-hand side IF (istage > 1) THEN DO j = 1, istage-1 ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) - CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z(1:N,j) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) IF (StartNewton) THEN - CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + rkAlpha(istage,j) * Z(1:N,j) END IF END DO END IF @@ -618,13 +616,13 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side - CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi + TMP(1:N) = Y(1:N) + Z(1:N,istage) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,DZ) ! DZ <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*DZ(1:N) - CALL WSCAL(N, H*rkGamma, DZ, 1) - CALL WAXPY (N, -ONE, Z(1,istage), 1, DZ, 1) - CALL WAXPY (N, ONE, G,1, DZ,1) + DZ(1:N) = DZ(1:N) * (H*rkGamma) + DZ(1:N) = DZ(1:N) - Z(1:N,istage) + DZ(1:N) = DZ(1:N) + G(1:N) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, DZ ) @@ -653,7 +651,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+DZ(:) - CALL WAXPY(N,ONE,DZ,1,Z(1,istage),1) + Z(1:N,istage) = Z(1:N,istage) + DZ(1:N) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) @@ -693,7 +691,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 - CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z_tlm(1:N,j,itlm) END DO END IF CALL SDIRK_Solve ( H, N, E_tlm, IP_tlm, IER, G ) @@ -712,7 +710,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) NewtonRate = MAX(NewtonRate,Roundoff)**0.8d0 !~~~> Starting values for Newton iterations - CALL Set2zero(N,Z_tlm(1,istage,itlm)) + Z_tlm(1:N,istage,itlm) = 0.0_dp !~~~> Prepare the loop-independent part of the right-hand side #ifdef FULL_ALGEBRA @@ -725,7 +723,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 - CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) + G(1:N) = G(1:N) + rkTheta(istage,j) * Z_tlm(1:N,j,itlm) END DO END IF @@ -778,7 +776,7 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) END IF !(TLMNewtonEst) ! Update solution: Z_tlm(:) <-- Z_tlm(:)+DZ(:) - CALL WAXPY(N,ONE,DZ,1,Z_tlm(1,istage,itlm),1) + Z_tlm(1:N,istage,itlm) = Z_tlm(1:N,istage,itlm) + DZ(1:N) ! Check error in Newton iterations IF (TLMNewtonEst) THEN @@ -809,19 +807,20 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 - CALL Set2zero(N,Yerr) + Yerr(1:N) = 0.0_dp DO i = 1,rkS - IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,Yerr,1) + IF (rkE(i)/=ZERO) Yerr(1:N) = Yerr(1:N) + rkE(i) * Z(1:N,i) END DO CALL SDIRK_Solve ( H, N, E, IP, IER, Yerr ) CALL SDIRK_ErrorNorm(N, Yerr, SCAL, Err) IF (TLMtruncErr) THEN - CALL Set2zero(NVAR*NTLM,Yerr_tlm) + Yerr_tlm(1:N,1:NTLM) = 0.0_dp DO itlm=1,NTLM DO j=1,rkS - IF (rkE(j) /= ZERO) CALL WAXPY(N,rkE(j),Z_tlm(1,j,itlm),1,Yerr_tlm(1,itlm),1) + IF (rkE(j) /= ZERO) Yerr_tlm(1:N,itlm) = Yerr_tlm(1:N,itlm) & + + rkE(j) * Z_tlm(1:N,j,itlm) END DO CALL SDIRK_Solve (H, N, E, IP, IER, Yerr_tlm(1,itlm)) END DO @@ -846,9 +845,10 @@ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) DO i = 1,rkS IF (rkD(i)/=ZERO) THEN - CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) + Y(1:N) = Y(1:N) + rkD(i) * Z(1:N,i) DO itlm = 1, NTLM - CALL WAXPY(N,rkD(i),Z_tlm(1,i,itlm),1,Y_tlm(1,itlm),1) + Y_tlm(1:N,itlm) = Y_tlm(1:N,itlm) & + + rkD(i) * Z_tlm(1:N,i,itlm) END DO END IF END DO @@ -1081,10 +1081,10 @@ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) KPP_REAL, INTENT(IN) :: E(LU_NONZERO) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) - KPP_REAL :: HGammaInv - HGammaInv = ONE/(H*rkGamma) - CALL WSCAL(N,HGammaInv,RHS,1) + ! This code replicates the output of the previous + ! call to WAXPY (@yantosca, 16 Oct 2025) + RHS(1:N) = RHS(1:N) * (ONE / (H*rkGamma)) #ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else diff --git a/int/seulex.f90 b/int/seulex.f90 index 2b0f4fe..348f090 100644 --- a/int/seulex.f90 +++ b/int/seulex.f90 @@ -426,7 +426,7 @@ SUBROUTINE ATMSEULEX( N,Tinitial,Tfinal,Y,H,RelTol,AbsTol, & !~~~> Unit roundoff (1+Roundoff>1) - Roundoff = WLAMCH('E') + Roundoff = EPSILON( 0.0_dp ) !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN diff --git a/src/gdata.h b/src/gdata.h index 8e970d8..36a2655 100644 --- a/src/gdata.h +++ b/src/gdata.h @@ -31,7 +31,7 @@ // Version numbers must be synchronized in CHANGELOG.md, src/gdata.h, // docs/source/conf.py and https://en.wikipedia.org/wiki/Kinetic_PreProcessor -#define KPP_VERSION "3.3.0" +#define KPP_VERSION "3.3.1" #ifndef _GDATA_H_ #define _GDATA_H_ diff --git a/util/Makefile_f90 b/util/Makefile_f90 index d4014fd..cbce72c 100644 --- a/util/Makefile_f90 +++ b/util/Makefile_f90 @@ -145,9 +145,9 @@ ifeq ($(HAS_STM_SP),1) STMOBJ += KPP_ROOT_StoichiomSP.o STMSPOBJ += KPP_ROOT_StoichiomSP.o endif -ifeq ($(HAS_HES),1) - HESSRC += KPP_ROOT_Stoichiom.f90 - HESOBJ += KPP_ROOT_Stoichiom.o +ifeq ($(HAS_STM),1) + STMSRC += KPP_ROOT_Stoichiom.f90 + STMOBJ += KPP_ROOT_Stoichiom.o endif UTLSRC = KPP_ROOT_Rates.f90 KPP_ROOT_Util.f90 KPP_ROOT_Monitor.f90 diff --git a/util/Makefile_upper_F90 b/util/Makefile_upper_F90 index fe6cc49..8780233 100644 --- a/util/Makefile_upper_F90 +++ b/util/Makefile_upper_F90 @@ -145,9 +145,9 @@ ifeq ($(HAS_STM_SP),1) STMOBJ += KPP_ROOT_StoichiomSP.o STMSPOBJ += KPP_ROOT_StoichiomSP.o endif -ifeq ($(HAS_HES),1) - HESSRC += KPP_ROOT_Stoichiom.F90 - HESOBJ += KPP_ROOT_Stoichiom.o +ifeq ($(HAS_STM),1) + STMSRC += KPP_ROOT_Stoichiom.F90 + STMOBJ += KPP_ROOT_Stoichiom.o endif UTLSRC = KPP_ROOT_Rates.F90 KPP_ROOT_Util.F90 KPP_ROOT_Monitor.F90 diff --git a/util/blas.f90 b/util/blas.f90 index cd8f337..c1a9a9b 100644 --- a/util/blas.f90 +++ b/util/blas.f90 @@ -8,320 +8,23 @@ ! Virginia Polytechnic Institute and State University !-------------------------------------------------------------- +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +!%%% NOTE: The following BLAS functions have been removed, %%% +!%%% as they now have been replaced by pure F90 code %%% +!%%% in the various integrator modules; %%% +!%%% %%% +!%%% (1) WCOPY %%% +!%%% (2) WAXPY %%% +!%%% (3) WSCAL %%% +!%%% (4) WLAMCH %%% +!%%% (5) WLAMCH_ADD %%% +!%%% (6) SET2ZERO %%% +!%%% (7) WADD %%% +!%%% (8) WDOT %%% +!%%% %%% +!%%% @yantosca, 17 Oct 2025 %%% +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!-------------------------------------------------------------- - SUBROUTINE WCOPY(N,X,incX,Y,incY) -!-------------------------------------------------------------- -! copies a vector, x, to a vector, y: y <- x -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1) -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - INTEGER :: i,incX,incY,M,MP1,N - KPP_REAL :: X(N),Y(N) - - IF (N.LE.0) RETURN - - M = MOD(N,8) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = X(i) - END DO - IF( N .LT. 8 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,8 - Y(i) = X(i) - Y(i + 1) = X(i + 1) - Y(i + 2) = X(i + 2) - Y(i + 3) = X(i + 3) - Y(i + 4) = X(i + 4) - Y(i + 5) = X(i + 5) - Y(i + 6) = X(i + 6) - Y(i + 7) = X(i + 7) - END DO - - END SUBROUTINE WCOPY - - -!-------------------------------------------------------------- - SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY) -!-------------------------------------------------------------- -! constant times a vector plus a vector: y <- y + Alpha*x -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1) -!-------------------------------------------------------------- - - INTEGER :: i,incX,incY,M,MP1,N - KPP_REAL :: X(N),Y(N),Alpha - KPP_REAL, PARAMETER :: ZERO = 0.0_dp - - IF (Alpha .EQ. ZERO) RETURN - IF (N .LE. 0) RETURN - - M = MOD(N,4) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = Y(i) + Alpha*X(i) - END DO - IF( N .LT. 4 ) RETURN - END IF - MP1 = M + 1 - DO i = MP1,N,4 - Y(i) = Y(i) + Alpha*X(i) - Y(i + 1) = Y(i + 1) + Alpha*X(i + 1) - Y(i + 2) = Y(i + 2) + Alpha*X(i + 2) - Y(i + 3) = Y(i + 3) + Alpha*X(i + 3) - END DO - - END SUBROUTINE WAXPY - - - -!-------------------------------------------------------------- - SUBROUTINE WSCAL(N,Alpha,X,incX) -!-------------------------------------------------------------- -! constant times a vector: x(1:N) <- Alpha*x(1:N) -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1) -!-------------------------------------------------------------- - - INTEGER :: i,incX,M,MP1,N - KPP_REAL :: X(N),Alpha - KPP_REAL, PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp - - IF (Alpha .EQ. ONE) RETURN - IF (N .LE. 0) RETURN - - M = MOD(N,5) - IF( M .NE. 0 ) THEN - IF (Alpha .EQ. (-ONE)) THEN - DO i = 1,M - X(i) = -X(i) - END DO - ELSEIF (Alpha .EQ. ZERO) THEN - DO i = 1,M - X(i) = ZERO - END DO - ELSE - DO i = 1,M - X(i) = Alpha*X(i) - END DO - END IF - IF( N .LT. 5 ) RETURN - END IF - MP1 = M + 1 - IF (Alpha .EQ. (-ONE)) THEN - DO i = MP1,N,5 - X(i) = -X(i) - X(i + 1) = -X(i + 1) - X(i + 2) = -X(i + 2) - X(i + 3) = -X(i + 3) - X(i + 4) = -X(i + 4) - END DO - ELSEIF (Alpha .EQ. ZERO) THEN - DO i = MP1,N,5 - X(i) = ZERO - X(i + 1) = ZERO - X(i + 2) = ZERO - X(i + 3) = ZERO - X(i + 4) = ZERO - END DO - ELSE - DO i = MP1,N,5 - X(i) = Alpha*X(i) - X(i + 1) = Alpha*X(i + 1) - X(i + 2) = Alpha*X(i + 2) - X(i + 3) = Alpha*X(i + 3) - X(i + 4) = Alpha*X(i + 4) - END DO - END IF - - END SUBROUTINE WSCAL - -!-------------------------------------------------------------- - KPP_REAL FUNCTION WLAMCH( C ) -!-------------------------------------------------------------- -! returns epsilon machine -! after LAPACK -! replace this by the function from the optimized LAPACK implementation: -! CALL SLAMCH('E') or CALL DLAMCH('E') -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - CHARACTER :: C - INTEGER :: i - KPP_REAL, SAVE :: Eps - KPP_REAL :: Suma - KPP_REAL, PARAMETER :: ONE=1.0_dp, HALF=0.5_dp - LOGICAL, SAVE :: First=.TRUE. - -!$OMP THREADPRIVATE( Eps, First ) - - IF (First) THEN - First = .FALSE. - Eps = HALF**(16) - DO i = 17, 80 - Eps = Eps*HALF - CALL WLAMCH_ADD(ONE,Eps,Suma) - IF (Suma.LE.ONE) GOTO 10 - END DO - PRINT*,'ERROR IN WLAMCH. EPS < ',Eps - RETURN -10 Eps = Eps*2 - i = i-1 - END IF - - WLAMCH = Eps - - END FUNCTION WLAMCH - - SUBROUTINE WLAMCH_ADD( A, B, Suma ) -! USE KPP_ROOT_Precision - - KPP_REAL A, B, Suma - Suma = A + B - - END SUBROUTINE WLAMCH_ADD -!-------------------------------------------------------------- - - -!-------------------------------------------------------------- - SUBROUTINE SET2ZERO(N,Y) -!-------------------------------------------------------------- -! copies zeros into the vector y: y <- 0 -! after BLAS -!-------------------------------------------------------------- - - INTEGER :: i,M,MP1,N - KPP_REAL :: Y(N) - KPP_REAL, PARAMETER :: ZERO = 0.0d0 - - IF (N.LE.0) RETURN - - M = MOD(N,8) - IF( M .NE. 0 ) THEN - DO i = 1,M - Y(i) = ZERO - END DO - IF( N .LT. 8 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,8 - Y(i) = ZERO - Y(i + 1) = ZERO - Y(i + 2) = ZERO - Y(i + 3) = ZERO - Y(i + 4) = ZERO - Y(i + 5) = ZERO - Y(i + 6) = ZERO - Y(i + 7) = ZERO - END DO - - END SUBROUTINE SET2ZERO - - -!-------------------------------------------------------------- - KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) -!-------------------------------------------------------------- -! dot produce: wdot = x(1:N)*y(1:N) -! only for incX=incY=1 -! after BLAS -! replace this by the function from the optimized BLAS implementation: -! CALL SDOT(N,X,1,Y,1) or CALL DDOT(N,X,1,Y,1) -!-------------------------------------------------------------- -! USE messy_mecca_kpp_Precision -!-------------------------------------------------------------- - IMPLICIT NONE - INTEGER :: N, incX, incY - KPP_REAL :: DX(N), DY(N) - - INTEGER :: i, IX, IY, M, MP1, NS - - WDOT = 0.0D0 - IF (N .LE. 0) RETURN - IF (incX .EQ. incY) IF (incX-1) 5,20,60 -! -! Code for unequal or nonpositive increments. -! - 5 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 - WDOT = WDOT + DX(IX)*DY(IY) - IX = IX + incX - IY = IY + incY - END DO - RETURN -! -! Code for both increments equal to 1. -! -! Clean-up loop so remaining vector length is a multiple of 5. -! - 20 M = MOD(N,5) - IF (M .EQ. 0) GO TO 40 - DO i = 1,M - WDOT = WDOT + DX(i)*DY(i) - END DO - IF (N .LT. 5) RETURN - 40 MP1 = M + 1 - DO i = MP1,N,5 - WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + & - DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4) - END DO - RETURN -! -! Code for equal, positive, non-unit increments. -! - 60 NS = N*incX - DO i = 1,NS,incX - WDOT = WDOT + DX(i)*DY(i) - END DO - - END FUNCTION WDOT - - -!-------------------------------------------------------------- - SUBROUTINE WADD(N,X,Y,Z) -!-------------------------------------------------------------- -! adds two vectors: z <- x + y -! BLAS - like -!-------------------------------------------------------------- -! USE KPP_ROOT_Precision - - INTEGER :: i, M, MP1, N - KPP_REAL :: X(N),Y(N),Z(N) - - IF (N.LE.0) RETURN - - M = MOD(N,5) - IF( M /= 0 ) THEN - DO i = 1,M - Z(i) = X(i) + Y(i) - END DO - IF( N < 5 ) RETURN - END IF - MP1 = M+1 - DO i = MP1,N,5 - Z(i) = X(i) + Y(i) - Z(i + 1) = X(i + 1) + Y(i + 1) - Z(i + 2) = X(i + 2) + Y(i + 2) - Z(i + 3) = X(i + 3) + Y(i + 3) - Z(i + 4) = X(i + 4) + Y(i + 4) - END DO - - END SUBROUTINE WADD - - - !-------------------------------------------------------------- SUBROUTINE WGEFA(N,A,Ipvt,info) !-------------------------------------------------------------- @@ -333,7 +36,7 @@ SUBROUTINE WGEFA(N,A,Ipvt,info) INTEGER :: N,Ipvt(N),info KPP_REAL :: A(N,N) KPP_REAL :: t, dmax, da - INTEGER :: j,k,l + INTEGER :: i,j,k,l KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0 info = 0 @@ -362,13 +65,19 @@ SUBROUTINE WGEFA(N,A,Ipvt,info) t = A(l,k); A(l,k) = A(k,k); A(k,k) = t END IF t = -ONE/A(k,k) - CALL WSCAL(n-k,t,A(k+1,k),1) +! CALL WSCAL(n-k,t,A(k+1,k),1) + DO i = k+1, n + A(i,k) = t * A(i,k) + END DO DO j = k+1, n t = A(l,j) IF (l /= k) THEN A(l,j) = A(k,j); A(k,j) = t END IF - CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) + !CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) + DO i = k+1, n + A(i,j) = A(i,j) + t * A(i,k) + END DO END DO END IF @@ -398,7 +107,7 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) CHARACTER :: trans KPP_REAL :: A(N,N),b(N) KPP_REAL :: t - INTEGER :: k,kb,l + INTEGER :: i, k,kb,l SELECT CASE (Trans) @@ -414,7 +123,10 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) b(l) = b(k) b(k) = t END IF - CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) + !CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) + DO i = k+1, n + b(i) = b(i) + t * a(i,k) + END DO END DO END IF ! now solve U*x = y @@ -422,21 +134,24 @@ SUBROUTINE WGESL(Trans,N,A,Ipvt,b) k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) - CALL WAXPY(k-1,t,a(1,k),1,b(1),1) + !CALL WAXPY(k-1,t,a(1,k),1,b(1),1) + DO i = 1, k-1 + b(i) = b(i) + t * a(i,k) + END DO END DO CASE ('t','T') ! Solve transpose(A) * x = b ! first solve trans(U)*y = b DO k = 1, n - t = WDOT(k-1,a(1,k),1,b(1),1) + t = DOT_PRODUCT( a(1:k-1, k), b(1:k-1) ) b(k) = (b(k) - t)/a(k,k) END DO ! now solve trans(L)*x = y IF (n >= 2) THEN DO kb = 1, n-1 k = n - kb - b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1) + b(k) = b(k) + DOT_PRODUCT( a(k+1:n, k), b(k+1:n) ) l = Ipvt(k) IF (l /= k) THEN t = b(l); b(l) = b(k); b(k) = t