diff --git a/.cmake-format.py b/.cmake-format.py new file mode 100644 index 00000000..6f2c0fee --- /dev/null +++ b/.cmake-format.py @@ -0,0 +1,243 @@ +# ---------------------------------- +# Options affecting listfile parsing +# ---------------------------------- +with section("parse"): + + # Specify structure for custom cmake functions + additional_commands = { + 'setup_test': { 'kwargs': {'GPUONLY': '0', 'LABELS': '+', 'MPI': '0', 'SOURCES': '+'}, + 'pargs': {'flags': ['2'], 'nargs': '0'}}, + } + + # Override configurations per-command where available + override_spec = {} + + # Specify variable tags. + vartags = [] + + # Specify property tags. + proptags = [] + +# ----------------------------- +# Options affecting formatting. +# ----------------------------- +with section("format"): + + # Disable formatting entirely, making cmake-format a no-op + disable = False + + # How wide to allow formatted cmake files + line_width = 100 + + # How many spaces to tab for indent + tab_size = 2 + + # If true, lines are indented using tab characters (utf-8 0x09) instead of + # space characters (utf-8 0x20). In cases where the layout would + # require a fractional tab character, the behavior of the fractional + # indentation is governed by + use_tabchars = False + + # If is True, then the value of this variable indicates how + # fractional indentions are handled during whitespace replacement. If set to + # 'use-space', fractional indentation is left as spaces (utf-8 0x20). If set + # to `round-up` fractional indentation is replaced with a single tab character + # (utf-8 0x09) effectively shifting the column to the next tabstop + fractional_tab_policy = 'use-space' + + # If an argument group contains more than this many sub-groups (parg or kwarg + # groups) then force it to a vertical layout. + max_subgroups_hwrap = 2 + + # If a positional argument group contains more than this many arguments, then + # force it to a vertical layout. + max_pargs_hwrap = 6 + + # If a cmdline positional group consumes more than this many lines without + # nesting, then invalidate the layout (and nest) + max_rows_cmdline = 2 + + # If true, separate flow control names from their parentheses with a space + separate_ctrl_name_with_space = False + + # If true, separate function names from parentheses with a space + separate_fn_name_with_space = False + + # If a statement is wrapped to more than one line, than dangle the closing + # parenthesis on its own line. + dangle_parens = True + + # If the trailing parenthesis must be 'dangled' on its on line, then align it + # to this reference: `prefix`: the start of the statement, `prefix-indent`: + # the start of the statement, plus one indentation level, `child`: align to + # the column of the arguments + dangle_align = 'prefix' + + # If the statement spelling length (including space and parenthesis) is + # smaller than this amount, then force reject nested layouts. + min_prefix_chars = 4 + + # If the statement spelling length (including space and parenthesis) is larger + # than the tab width by more than this amount, then force reject un-nested + # layouts. + max_prefix_chars = 10 + + # If a candidate layout is wrapped horizontally but it exceeds this many + # lines, then reject the layout. + max_lines_hwrap = 2 + + # What style line endings to use in the output. + line_ending = 'unix' + + # Format command names consistently as 'lower' or 'upper' case + command_case = 'canonical' + + # Format keywords consistently as 'lower' or 'upper' case + keyword_case = 'unchanged' + + # A list of command names which should always be wrapped + always_wrap = [] + + # If true, the argument lists which are known to be sortable will be sorted + # lexicographicall + enable_sort = True + + # If true, the parsers may infer whether or not an argument list is sortable + # (without annotation). + autosort = False + + # By default, if cmake-format cannot successfully fit everything into the + # desired linewidth it will apply the last, most aggressive attempt that it + # made. If this flag is True, however, cmake-format will print error, exit + # with non-zero status code, and write-out nothing + require_valid_layout = False + + # A dictionary mapping layout nodes to a list of wrap decisions. See the + # documentation for more information. + layout_passes = {} + +# ------------------------------------------------ +# Options affecting comment reflow and formatting. +# ------------------------------------------------ +with section("markup"): + + # What character to use for bulleted lists + bullet_char = '*' + + # What character to use as punctuation after numerals in an enumerated list + enum_char = '.' + + # If comment markup is enabled, don't reflow the first comment block in each + # listfile. Use this to preserve formatting of your copyright/license + # statements. + first_comment_is_literal = False + + # If comment markup is enabled, don't reflow any comment block which matches + # this (regex) pattern. Default is `None` (disabled). + literal_comment_pattern = None + + # Regular expression to match preformat fences in comments default= + # ``r'^\s*([`~]{3}[`~]*)(.*)$'`` + fence_pattern = '^\\s*([`~]{3}[`~]*)(.*)$' + + # Regular expression to match rulers in comments default= + # ``r'^\s*[^\w\s]{3}.*[^\w\s]{3}$'`` + ruler_pattern = '^\\s*[^\\w\\s]{3}.*[^\\w\\s]{3}$' + + # If a comment line matches starts with this pattern then it is explicitly a + # trailing comment for the preceding argument. Default is '#<' + explicit_trailing_pattern = '#<' + + # If a comment line starts with at least this many consecutive hash + # characters, then don't lstrip() them off. This allows for lazy hash rulers + # where the first hash char is not separated by space + hashruler_min_length = 10 + + # If true, then insert a space between the first hash char and remaining hash + # chars in a hash ruler, and normalize its length to fill the column + canonicalize_hashrulers = True + + # enable comment markup parsing and reflow + enable_markup = True + +# ---------------------------- +# Options affecting the linter +# ---------------------------- +with section("lint"): + + # a list of lint codes to disable + disabled_codes = [] + + # regular expression pattern describing valid function names + function_pattern = '[0-9a-z_]+' + + # regular expression pattern describing valid macro names + macro_pattern = '[0-9A-Z_]+' + + # regular expression pattern describing valid names for variables with global + # (cache) scope + global_var_pattern = '[A-Z][0-9A-Z_]+' + + # regular expression pattern describing valid names for variables with global + # scope (but internal semantic) + internal_var_pattern = '_[A-Z][0-9A-Z_]+' + + # regular expression pattern describing valid names for variables with local + # scope + local_var_pattern = '[a-z][a-z0-9_]+' + + # regular expression pattern describing valid names for privatedirectory + # variables + private_var_pattern = '_[0-9a-z_]+' + + # regular expression pattern describing valid names for public directory + # variables + public_var_pattern = '[A-Z][0-9A-Z_]+' + + # regular expression pattern describing valid names for function/macro + # arguments and loop variables. + argument_var_pattern = '[a-z][a-z0-9_]+' + + # regular expression pattern describing valid names for keywords used in + # functions or macros + keyword_pattern = '[A-Z][0-9A-Z_]+' + + # In the heuristic for C0201, how many conditionals to match within a loop in + # before considering the loop a parser. + max_conditionals_custom_parser = 2 + + # Require at least this many newlines between statements + min_statement_spacing = 1 + + # Require no more than this many newlines between statements + max_statement_spacing = 2 + max_returns = 6 + max_branches = 12 + max_arguments = 5 + max_localvars = 15 + max_statements = 50 + +# ------------------------------- +# Options affecting file encoding +# ------------------------------- +with section("encode"): + + # If true, emit the unicode byte-order mark (BOM) at the start of the file + emit_byteorder_mark = False + + # Specify the encoding of the input file. Defaults to utf-8 + input_encoding = 'utf-8' + + # Specify the encoding of the output file. Defaults to utf-8. Note that cmake + # only claims to support utf-8 so be careful when using anything else + output_encoding = 'utf-8' + +# ------------------------------------- +# Miscellaneous configurations options. +# ------------------------------------- +with section("misc"): + + # A dictionary containing any per-command configuration overrides. Currently + # only `command_case` is supported. + per_command = {} + diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9aa7e602..e4b94dd3 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -19,6 +19,18 @@ MALLOC_PERTURB_: 1 +check-format: + stage: .pre + image: $CI_REGISTRY/nexgf/libnegf/check-format + script: + - >- + find . -name CMakeLists.txt -print0 + | xargs -0 -L1 -- cmake-format --check -- + - >- + find . -iname '*.cu' -o -iname '*.cpp' -print0 + | xargs -0 -L1 -- clang-format --dry-run -Werror -- + + build: extends: .setup stage: build @@ -47,4 +59,5 @@ test: stage: test script: - cd "$BUILD_DIR" - - ctest -L normal --output-on-failure + - ctest --output-on-failure -LE long + - ctest --output-on-failure -L long diff --git a/CMakeLists.txt b/CMakeLists.txt index da5accb0..cc04720b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,9 @@ cmake_minimum_required(VERSION 3.18) +if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.31) + cmake_policy(SET CMP0177 NEW) +endif() + list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/cmake) include(LibNegfUtils) libnegf_load_build_settings() @@ -7,37 +11,41 @@ libnegf_load_build_settings() # Some setting output message(STATUS "libNEGF WITH_MPI: " ${WITH_MPI}) message(STATUS "libNEGF WITH_GPU: " ${WITH_TRANSPORT_GPU}) -#message(STATUS "libNEGF WITH_HILBERT: " ${WITH_HILBERT}) +# message(STATUS "libNEGF WITH_HILBERT: " ${WITH_HILBERT}) if(WITH_TRANSPORT_GPU) - project(libNEGF VERSION 1.2.1 LANGUAGES Fortran C CUDA CXX) + project( + libNEGF + VERSION 1.2.1 + LANGUAGES Fortran C CUDA CXX + ) else() - project(libNEGF VERSION 1.2.1 LANGUAGES Fortran C CXX) + project( + libNEGF + VERSION 1.2.1 + LANGUAGES Fortran C CXX + ) endif() - # disable compiler-specific extensions set(CMAKE_C_STANDARD 99) set(CMAKE_C_STANDARD_REQUIRED ON) -# Note: enforcing a standard in Fortran compilers has a different effect that in C compilers, c.f. [cmake#22235](https://gitlab.kitware.com/cmake/cmake/-/issues/22235). - +# Note: enforcing a standard in Fortran compilers has a different effect that in C compilers, c.f. +# [cmake#22235](https://gitlab.kitware.com/cmake/cmake/-/issues/22235). include(externalMpifx) include(CMakePackageConfigHelpers) - # find a Python3 interpreter for FYPP include(FindPython3) find_package(Python3 REQUIRED COMPONENTS Interpreter) - set(CMAKE_PREFIX_PATH "${CMAKE_CURRENT_SOURCE_DIR}/ext_fypp/;${CMAKE_PREFIX_PATH}") find_program(FYPP "fypp") include(GNUInstallDirs) - if(LAPACK_LIBRARIES AND LAPACK_LIBRARY_DIRS) message(STATUS "libNEGF LAPACK DIR: " ${LAPACK_LIBRARY_DIRS}) message(STATUS "libNEGF BLAS: " ${BLAS_LIBRARY}) @@ -53,7 +61,6 @@ if(WITH_MPI) find_or_fetch_mpifx() endif() - # Subdirectories. add_subdirectory(ext_system) add_subdirectory(src) @@ -64,12 +71,13 @@ if(BUILD_TESTING) add_subdirectory(tests) endif() - # # Installation # -set(CMAKE_INSTALL_PREFIX "/usr/local/libnegf${libNEGF_VERSION_MAJOR}.${libNEGF_VERSION_MINOR}" CACHE PATH - "default installation path") +set(CMAKE_INSTALL_PREFIX + "/usr/local/libnegf${libNEGF_VERSION_MAJOR}.${libNEGF_VERSION_MINOR}" + CACHE PATH "default installation path" +) message(STATUS "Default installation path: " ${CMAKE_INSTALL_PREFIX}) @@ -77,23 +85,27 @@ add_library(Negf INTERFACE) target_link_libraries(Negf INTERFACE negf) install(TARGETS Negf EXPORT negf-targets) -install(EXPORT negf-targets +install( + EXPORT negf-targets FILE negf-targets.cmake NAMESPACE Negf:: DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/negf" - EXPORT_LINK_INTERFACE_LIBRARIES) + EXPORT_LINK_INTERFACE_LIBRARIES +) configure_package_config_file( ${CMAKE_CURRENT_SOURCE_DIR}/utils/export/negf-config.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config.cmake - INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/negf) + INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/negf +) write_basic_package_version_file( ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config-version.cmake VERSION ${PROJECT_VERSION} - COMPATIBILITY SameMajorVersion) + COMPATIBILITY SameMajorVersion +) -install( - FILES ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config.cmake - ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config-version.cmake - DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/negf) +install(FILES ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config.cmake + ${CMAKE_CURRENT_BINARY_DIR}/cmake/negf-config-version.cmake + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/negf +) diff --git a/README.adoc b/README.adoc index 49369adb..90f47549 100644 --- a/README.adoc +++ b/README.adoc @@ -1,5 +1,4 @@ -libNEGF -======= += libNEGF libNEGF is a general library for Non Equilibrium Green's Functions. @@ -37,6 +36,86 @@ Running the tests additionally requires git-lfs for downloading test inputs. The extension of test input files is `.dat`. +== Running the Tests on a Supercomputer + +On some HPC systems, front-end nodes do not feature GPUs (e.g., on LUMI) or the +front-end nodes may have an instruction set architecture different from compute +nodes (e.g., on Fugaku with its Intel front-end and Arm64FX compute nodes). On +these systems, it is best to submit all tests to the batch scheduler. + +This can be achieved with the following steps. First, determine the account name and a queue (or _partition_) for job submission. Next, select a time limit. Tests that do not possess the label `long` run quickly even on personal computers. Therefore, we suggest a time limit of at most five minutes (this is per task). Finally, determine the number of MPI tasks within the job. Our recommendation is to run on at most one node with one job per NUMA domain. + +CAUTION: Avoid large numbers of MPI tasks (e.g., by launching one task per virtual CPU core). The tests may not scale very well. + +With the values above in mind, enable `LIBNEGF_TEST_WITH_MPIEXEC`, set +`MPIEXEC_EXECUTABLE` to the absolute path to the batch scheduler executable, and +have `MPIEXEC_PREFLAGS` contain all the batch scheduler arguments as a CMake +list; in a CMake list, list items are separated by a semicolon. Here is an +example for Slurm: + +[source,shell,linenums] +---- +cmake \ + -DLIBNEGF_TEST_WITH_MPIEXEC=ON \ + -DMPIEXEC_EXECUTABLE="$(which srun)" \ + -DMPIEXEC_NUMPROC_FLAG='--ntasks' \ + -DMPIEXEC_MAX_NUMPROCS=4 \ + -DMPIEXEC_PREFLAGS='--account=mat4energy;--partition=develbooster;--nodes=1;--gpus-per-task=1;--time=1' \ + ... +---- +For faster execution, the tests can be run in an existing job allocation, e.g., +by calling `salloc` (for Slurm) or by using a reservation. Within the allocation, just call `ctest`. + +Note the following: + +* A time limit of one minute is sufficient for tests without `long` tag on JUWELS Booster and LUMI. +* The time limit of the job and the time limit of test CMake test submitting the job are two different things. +* Redirecting the output to files (e.g., with Slurm with `--output=stdout.txt`) breaks the mpifx tests because these tests examine the standard output instead of the exit status. + + + +== Contributing + +libNEGF enforces a coding style in C, C++, and CMake. Ideally, all committed code should contain only code adhering to the style guidelines. To support contributors in this regard, a sample git pre-commit hook checking the format is provided in link:git-pre-commit-hook.sh[`git-pre-commit-hook.sh`]. To install the file, append its contents to `.git/hooks/pre-commit` and make the latter file executable; the script requires cmake-format and clang-format to be installed. See the https://git-scm.com/book/en/v2/Customizing-Git-Git-Hooks[Git Book: _Customizing Git - Git Hooks_] or https://man7.org/linux/man-pages/man5/githooks.5.html[`man 5 githooks`] for more detailed instructions. The pre-commit script requires clang-format and cmake-format to be installed. + +[NOTE] +==== +Given the same source code and configuration, clang-format does _not_ promise to create the same formatting across different versions. The following clang-format major releases were tested in libnegf#117 and generate the same output unless noted otherwise: + +* 11 (Debian 11, *not* working) +* 14 (Debian 12) +* 18 (JUWELS Booster, Stages/2025) +* 19 (Debian 13) +* 20 (JUWELS Booster, Stages/2026) +==== + + + + +== Memory Scaling + +The peak memory consumption is a function of the device size and the number of +its principal layers (PLs) with all PLs assumed to be identical. This allows +one to estimate maximum problem sizes for a given set of resources. + +The principal layer size is +``` +PL_size = (N_atoms × N_orbitals)² × 16 bytes +``` +in double precision, where the number of atoms `N_atoms` scales quadratically +with the device size. The number of orbitals `N_orbitals` is none and includes +the d orbitals. `N_pl` denotes the number of principal layers. + +.libNEGF peak memory consumption +[options="header"] +|=== +| Computation | Peak memory consumption +| Coherent transmission, `computeLDOS=No` | 7 × PL_size +| Coherent transmission, `computeLDOS=Yes` | 3 × N_pl × PL_size +| Density matrix (coherent transmission) | 7 × N_pl × PL_size +| Inelastic code | 12 × N_pl × PL_size +|=== + == Generating libNEGF Input diff --git a/config.cmake b/config.cmake index b460c0c7..18ef7bf9 100644 --- a/config.cmake +++ b/config.cmake @@ -39,6 +39,22 @@ set(INSTALL_INCLUDE_DIR "" CACHE PATH set(INSTALL_MOD_DIR "${INSTALL_INCLUDE_DIR}/modfiles" CACHE PATH "Installation directory for Fortran module files (within standard include folder)") +# Parallel tests are launched with the command: +# env OMP_NUM_THREADS=${TEST_NUM_THREADS} +# ${MPIEXEC_EXECUTABLE} +# ${MPIEXEC_NUMPROC_FLAG} +# ${TEST_NUM_PROCS} ${MPIEXEC_PREFLAGS} +# ./${testname} ${MPIEXEC_POSTFLAGS} +# +# To run ctest in a HPC cluster set: +#set(MPIEXEC_EXECUTABLE "srun" CACHE STRING "Define the executable") +#set(MPIEXEC_NUMPROC_FLAG "--ntasks" CACHE STRING "Define parallel flag") +#set(MPIEXEC_PREFLAGS "" CACHE STRING "Define prepend flags") +#set(MPIEXEC_POSTFLAGS "" CACHE STRING "Define postpend flags") + +set(TEST_NUM_PROCS ${MPIEXEC_MAX_NUMPROCS} CACHE STRING "Define number of tasks for tests") +set(TEST_NUM_THREADS "1" CACHE STRING "Define OMP_NUM_THREADS for tests") + ####################################################################################### # Fortran and C compilers, manual setting example ####################################################################################### diff --git a/ext_system/CMakeLists.txt b/ext_system/CMakeLists.txt index adda84a1..660637a1 100644 --- a/ext_system/CMakeLists.txt +++ b/ext_system/CMakeLists.txt @@ -1,6 +1,4 @@ -set(sources - sys_calls.c - system_calls.f90) +set(sources sys_calls.c system_calls.f90) add_library(syscalls_objlib OBJECT ${sources}) @@ -9,9 +7,10 @@ set(BUILD_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/include) set_target_properties(syscalls_objlib PROPERTIES Fortran_MODULE_DIRECTORY ${BUILD_MOD_DIR}) set_target_properties(syscalls_objlib PROPERTIES POSITION_INDEPENDENT_CODE ${BUILD_SHARED_LIBS}) -target_include_directories(syscalls_objlib PUBLIC - $ - $) +target_include_directories( + syscalls_objlib PUBLIC $ + $ +) if(INSTALL_INCLUDE_FILES) install(DIRECTORY ${BUILD_MOD_DIR}/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${INSTALL_MOD_DIR}) diff --git a/git-pre-commit-hook.sh b/git-pre-commit-hook.sh new file mode 100644 index 00000000..f47e90b9 --- /dev/null +++ b/git-pre-commit-hook.sh @@ -0,0 +1,54 @@ +#!/bin/bash +# Author: Christoph Conrads (Forschungszentrum Jülich) +# +# The git pre-commit hook checks if Python and C++ code is properly formatted. +# +# The code requires YAPF for Python and clang-format for C/C++/CUDA. + +set -e +set -u + + +# check for presence of formatting tools +if ! which clang-format >/dev/null +then + >/dev/stderr echo 'clang-format not found. aborting.' + exit 1 +fi + +if ! python3 -c 'import yapf' 2>/dev/null +then + >/dev/stderr echo 'Python3 module YAPF not found. aborting.' + exit 1 +fi + + +needs_format=false + + +# check files +files="$(git diff --cached --name-only --diff-filter=ACM)" + +while IFS='' read -r file +do + if [[ "$file" =~ [.]([ch]pp|cu)$ ]] + then + if ! $(clang-format --dry-run --Werror -- "$file" 2>/dev/null) + then + echo "$file must be formatted with clang-format" + needs_format=true + fi + elif [[ "$file" =~ [.]py$ ]] + then + if ! $(python3 -m yapf --quiet -- "$file") + then + echo "$file must be formatted with YAPF" + needs_format=true + fi + fi +done <<<"$files" + +if [ "$needs_format" = true ] +then + exit 1 +fi diff --git a/include/libnegf/cublas_v2.h b/include/libnegf/cublas_v2.h index a805259e..6029e733 100644 --- a/include/libnegf/cublas_v2.h +++ b/include/libnegf/cublas_v2.h @@ -32,7 +32,7 @@ namespace libnegf { cublasStatus_t cublasAsum( - cublasHandle_t handle, size_t n, const cuComplex* x, size_t incx, + cublasHandle_t handle, int32_t n, const cuComplex* x, int32_t incx, float* result ) { assert(n <= INT_MAX); @@ -41,7 +41,7 @@ cublasStatus_t cublasAsum( } cublasStatus_t cublasAsum( - cublasHandle_t handle, size_t n, const cuDoubleComplex* x, size_t incx, + cublasHandle_t handle, int32_t n, const cuDoubleComplex* x, int32_t incx, double* result ) { assert(n <= INT_MAX); @@ -51,8 +51,8 @@ cublasStatus_t cublasAsum( cublasStatus_t cublasAxpy( - cublasHandle_t handle, size_t n, const cuComplex* alpha, const cuComplex* x, - size_t incx, cuComplex* y, size_t incy + cublasHandle_t handle, int32_t n, const cuComplex* alpha, const cuComplex* x, + int32_t incx, cuComplex* y, int32_t incy ) { assert(n <= INT_MAX); assert(alpha); @@ -63,8 +63,8 @@ cublasStatus_t cublasAxpy( } cublasStatus_t cublasAxpy( - cublasHandle_t handle, size_t n, const cuDoubleComplex* alpha, - const cuDoubleComplex* x, size_t incx, cuDoubleComplex* y, size_t incy + cublasHandle_t handle, int32_t n, const cuDoubleComplex* alpha, + const cuDoubleComplex* x, int32_t incx, cuDoubleComplex* y, int32_t incy ) { assert(n <= INT_MAX); assert(alpha); @@ -76,8 +76,8 @@ cublasStatus_t cublasAxpy( cublasStatus_t cublasCopy( - cublasHandle_t handle, size_t n, const cuComplex* d_x, size_t incx, - cuComplex* d_y, size_t incy + cublasHandle_t handle, int32_t n, const cuComplex* d_x, int32_t incx, + cuComplex* d_y, int32_t incy ) { assert(n <= INT_MAX); assert(d_x); @@ -87,8 +87,8 @@ cublasStatus_t cublasCopy( } cublasStatus_t cublasCopy( - cublasHandle_t handle, size_t n, const cuDoubleComplex* d_x, size_t incx, - cuDoubleComplex* d_y, size_t incy + cublasHandle_t handle, int32_t n, const cuDoubleComplex* d_x, int32_t incx, + cuDoubleComplex* d_y, int32_t incy ) { assert(n <= INT_MAX); assert(d_x); @@ -101,8 +101,8 @@ cublasStatus_t cublasCopy( // cublasXdot cublasStatus_t cublasDot( - cublasHandle_t handle, size_t n, const float* x, size_t incx, - const float* y, size_t incy, float* result + cublasHandle_t handle, int32_t n, const float* x, int32_t incx, + const float* y, int32_t incy, float* result ) { assert(n); assert(x); @@ -113,8 +113,8 @@ cublasStatus_t cublasDot( } cublasStatus_t cublasDot( - cublasHandle_t handle, size_t n, const double* x, size_t incx, - const double* y, size_t incy, double* result + cublasHandle_t handle, int32_t n, const double* x, int32_t incx, + const double* y, int32_t incy, double* result ) { assert(n); assert(x); @@ -129,8 +129,8 @@ cublasStatus_t cublasDot( cublasStatus_t cublasGeam( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, const float* alpha, const float* A, size_t lda, - const float* beta, const float* B, size_t ldb, float* C, size_t ldc + int32_t m, int32_t n, const float* alpha, const float* A, int32_t lda, + const float* beta, const float* B, int32_t ldb, float* C, int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -144,8 +144,8 @@ cublasStatus_t cublasGeam( cublasStatus_t cublasGeam( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, const double* alpha, const double* A, size_t lda, - const double* beta, const double* B, size_t ldb, double* C, size_t ldc + int32_t m, int32_t n, const double* alpha, const double* A, int32_t lda, + const double* beta, const double* B, int32_t ldb, double* C, int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -159,9 +159,9 @@ cublasStatus_t cublasGeam( cublasStatus_t cublasGeam( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, const cuComplex* alpha, const cuComplex* A, size_t lda, - const cuComplex* beta, const cuComplex* B, size_t ldb, cuComplex* C, - size_t ldc + int32_t m, int32_t n, const cuComplex* alpha, const cuComplex* A, int32_t lda, + const cuComplex* beta, const cuComplex* B, int32_t ldb, cuComplex* C, + int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -175,9 +175,9 @@ cublasStatus_t cublasGeam( cublasStatus_t cublasGeam( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, const cuDoubleComplex* alpha, const cuDoubleComplex* A, - size_t lda, const cuDoubleComplex* beta, const cuDoubleComplex* B, - size_t ldb, cuDoubleComplex* C, size_t ldc + int32_t m, int32_t n, const cuDoubleComplex* alpha, const cuDoubleComplex* A, + int32_t lda, const cuDoubleComplex* beta, const cuDoubleComplex* B, + int32_t ldb, cuDoubleComplex* C, int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -194,9 +194,9 @@ cublasStatus_t cublasGeam( cublasStatus_t cublasGemm( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, size_t k, const float* alpha, const float* A, - size_t lda, const float* B, size_t ldb, const float* beta, float* C, - size_t ldc + int32_t m, int32_t n, int32_t k, const float* alpha, const float* A, + int32_t lda, const float* B, int32_t ldb, const float* beta, float* C, + int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -211,9 +211,9 @@ cublasStatus_t cublasGemm( cublasStatus_t cublasGemm( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, size_t k, const double* alpha, const double* A, - size_t lda, const double* B, size_t ldb, const double* beta, double* C, - size_t ldc + int32_t m, int32_t n, int32_t k, const double* alpha, const double* A, + int32_t lda, const double* B, int32_t ldb, const double* beta, double* C, + int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -228,9 +228,9 @@ cublasStatus_t cublasGemm( cublasStatus_t cublasGemm( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, size_t k, const cuComplex* alpha, const cuComplex* A, - size_t lda, const cuComplex* B, size_t ldb, const cuComplex* beta, - cuComplex* C, size_t ldc + int32_t m, int32_t n, int32_t k, const cuComplex* alpha, const cuComplex* A, + int32_t lda, const cuComplex* B, int32_t ldb, const cuComplex* beta, + cuComplex* C, int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -245,9 +245,9 @@ cublasStatus_t cublasGemm( cublasStatus_t cublasGemm( cublasHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, - size_t m, size_t n, size_t k, const cuDoubleComplex* alpha, - const cuDoubleComplex* A, size_t lda, const cuDoubleComplex* B, size_t ldb, - const cuDoubleComplex* beta, cuDoubleComplex* C, size_t ldc + int32_t m, int32_t n, int32_t k, const cuDoubleComplex* alpha, + const cuDoubleComplex* A, int32_t lda, const cuDoubleComplex* B, int32_t ldb, + const cuDoubleComplex* beta, cuDoubleComplex* C, int32_t ldc ) { assert(m <= INT_MAX); assert(n <= INT_MAX); diff --git a/include/libnegf/cusolverDn.h b/include/libnegf/cusolverDn.h index 88fb8b8e..73c8d928 100644 --- a/include/libnegf/cusolverDn.h +++ b/include/libnegf/cusolverDn.h @@ -31,7 +31,7 @@ namespace libnegf { // the deprecated cuSolver function of the same name. cusolverStatus_t cusolverDngetrf( - cusolverDnHandle_t handle, size_t m, size_t n, cuComplex* A, size_t lda, cuComplex* workspace, int *d_ipiv, int* d_info) + cusolverDnHandle_t handle, int32_t m, int32_t n, cuComplex* A, int32_t lda, cuComplex* workspace, int *d_ipiv, int* d_info) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -45,7 +45,7 @@ cusolverStatus_t cusolverDngetrf( } cusolverStatus_t cusolverDngetrf( - cusolverDnHandle_t handle, size_t m, size_t n, cuDoubleComplex* A, size_t lda, cuDoubleComplex* workspace, int *d_ipiv, int* d_info) + cusolverDnHandle_t handle, int32_t m, int32_t n, cuDoubleComplex* A, int32_t lda, cuDoubleComplex* workspace, int *d_ipiv, int* d_info) { assert(m <= INT_MAX); assert(n <= INT_MAX); @@ -60,7 +60,7 @@ cusolverStatus_t cusolverDngetrf( cusolverStatus_t cusolverDngetrs( - cusolverDnHandle_t handle, cublasOperation_t trans, size_t n, size_t nrhs, const cuComplex* A, size_t lda, const int* d_ipiv, cuComplex* B, size_t ldb, int* d_info) + cusolverDnHandle_t handle, cublasOperation_t trans, int32_t n, int32_t nrhs, const cuComplex* A, int32_t lda, const int* d_ipiv, cuComplex* B, int32_t ldb, int* d_info) { assert(n <= INT_MAX); assert(nrhs <= INT_MAX); @@ -73,7 +73,7 @@ cusolverStatus_t cusolverDngetrs( } cusolverStatus_t cusolverDngetrs( - cusolverDnHandle_t handle, cublasOperation_t trans, size_t n, size_t nrhs, const cuDoubleComplex* A, size_t lda, const int* d_ipiv, cuDoubleComplex* B, size_t ldb, int* d_info) + cusolverDnHandle_t handle, cublasOperation_t trans, int32_t n, int32_t nrhs, const cuDoubleComplex* A, int32_t lda, const int* d_ipiv, cuDoubleComplex* B, int32_t ldb, int* d_info) { assert(n <= INT_MAX); assert(nrhs <= INT_MAX); @@ -87,7 +87,7 @@ cusolverStatus_t cusolverDngetrs( cusolverStatus_t cusolverDngetrf_bufferSize( - cusolverDnHandle_t handle, size_t m, size_t n, cuComplex* A, size_t lda, + cusolverDnHandle_t handle, int32_t m, int32_t n, cuComplex* A, int32_t lda, int* lwork ) { assert(m <= INT_MAX); @@ -100,8 +100,8 @@ cusolverStatus_t cusolverDngetrf_bufferSize( } cusolverStatus_t cusolverDngetrf_bufferSize( - cusolverDnHandle_t handle, size_t m, size_t n, cuDoubleComplex* A, - size_t lda, int* lwork + cusolverDnHandle_t handle, int32_t m, int32_t n, cuDoubleComplex* A, + int32_t lda, int* lwork ) { assert(m <= INT_MAX); assert(n <= INT_MAX); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6200f754..5898155a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,43 +1,44 @@ set(sources-fpp - clock.F90 - ln_extract.F90 - input_output.F90 - ln_cache.F90 - lib_param.F90 - ln_enums.F90 - ln_precision.F90 - mpi_globals.F90 - ln_messages.F90 - contselfenergy.F90 - inversions.F90 - libnegf.F90 - ln_structure.F90 - outmatrix.F90 - distributions.F90 - ln_allocation.F90 - load.F90 - sparsekit_drv.F90 - integrations.F90 - globals.F90 - iterative.F90 - ln_constants.F90 - mat_def.F90 - population.F90 - energy_mesh.F90 - equiv_kpoints.F90 - interactions.F90 - ln_elastic.F90 - ln_inelastic.F90 - elphdd.F90 - elphdb.F90 - elphds.F90 - elphinel.F90 - self_energy.F90 - scba.F90 - sparskit/skit_blassm.F90 - sparskit/skit_formats.F90 - sparskit/skit_module.F90 - sparskit/skit_unary.F90) + clock.F90 + ln_extract.F90 + input_output.F90 + ln_cache.F90 + lib_param.F90 + ln_enums.F90 + ln_precision.F90 + mpi_globals.F90 + ln_messages.F90 + contselfenergy.F90 + inversions.F90 + libnegf.F90 + ln_structure.F90 + outmatrix.F90 + distributions.F90 + ln_allocation.F90 + load.F90 + sparsekit_drv.F90 + integrations.F90 + globals.F90 + iterative.F90 + ln_constants.F90 + mat_def.F90 + population.F90 + energy_mesh.F90 + equiv_kpoints.F90 + interactions.F90 + ln_elastic.F90 + ln_inelastic.F90 + elphdd.F90 + elphdb.F90 + elphds.F90 + elphinel.F90 + self_energy.F90 + scba.F90 + sparskit/skit_blassm.F90 + sparskit/skit_formats.F90 + sparskit/skit_module.F90 + sparskit/skit_unary.F90 +) if(WITH_HILBERT) list(APPEND sources-fpp transform.F90) @@ -48,17 +49,15 @@ if(WITH_TRANSPORT_GPU) list(APPEND sources-fpp iterative_gpu.F90) set(sources-cuda gpuroutines.cu) set_source_files_properties( - gpuroutines.cu - PROPERTIES INCLUDE_DIRECTORIES ${PROJECT_SOURCE_DIR}/include) + gpuroutines.cu PROPERTIES INCLUDE_DIRECTORIES ${PROJECT_SOURCE_DIR}/include + ) else() list(APPEND sources-fpp iterative_cpu.F90) endif() -set(sources-f90 - api/libnegfAPICommon.f90 - api/libnegf_api.f90) +set(sources-f90 api/libnegfAPICommon.f90 api/libnegf_api.f90) -#execute_process(COMMAND git describe OUTPUT_VARIABLE gitrevision) +# execute_process(COMMAND git describe OUTPUT_VARIABLE gitrevision) set(gitrevision "000") string(TIMESTAMP compdate "%Y-%m-%d") @@ -67,7 +66,6 @@ if(WITH_MPI) list(APPEND fyppdefs -DMPI) endif() - if(WITH_TRANSPORT_GPU) list(APPEND fyppdefs -DGPU) endif() @@ -77,26 +75,25 @@ foreach(fppsrc IN LISTS sources-fpp) string(REGEX REPLACE "\\.F90" ".f90" f90src ${fppsrc}) add_custom_command( OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${f90src} - COMMAND ${FYPP} - -I${PROJECT_SOURCE_DIR}/include - ${fyppdefs} ${CMAKE_CURRENT_SOURCE_DIR}/${fppsrc} - ${CMAKE_CURRENT_BINARY_DIR}/${f90src} + COMMAND ${FYPP} -I${PROJECT_SOURCE_DIR}/include ${fyppdefs} + ${CMAKE_CURRENT_SOURCE_DIR}/${fppsrc} ${CMAKE_CURRENT_BINARY_DIR}/${f90src} MAIN_DEPENDENCY ${CMAKE_CURRENT_SOURCE_DIR}/${fppsrc} DEPENDS ${PROJECT_SOURCE_DIR}/include/assert.fypp - VERBATIM) + VERBATIM + ) list(APPEND sources-fpp-f90 ${CMAKE_CURRENT_BINARY_DIR}/${f90src}) endforeach() - - if(WITH_TRANSPORT_GPU) - add_library(negf ${sources-fpp-f90} ${sources-f90} ${sources-cuda}) - else() - add_library(negf ${sources-fpp-f90} ${sources-f90}) - endif() +if(WITH_TRANSPORT_GPU) + add_library(negf ${sources-fpp-f90} ${sources-f90} ${sources-cuda}) +else() + add_library(negf ${sources-fpp-f90} ${sources-f90}) +endif() target_sources(negf PRIVATE $) -target_include_directories(negf PUBLIC - $>) +target_include_directories( + negf PUBLIC $> +) set(BUILD_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/include) @@ -111,16 +108,13 @@ endif() if(WITH_MPI) find_package(MPI REQUIRED) target_link_libraries(negf PUBLIC MpiFx::MpiFx) - # MPI::MPI_Fortran cannot be made PUBLIC because the compile options of this - # target may cause errors or warnings when passed to a C compiler (e.g., when - # using MPICH 4.0.2 on Debian 12). Unfortunately we cannot manually add - # INTERFACE compile options because this leads to build problems on - # * JUWELS with Rocky Linux 8, ParaStationMPI 5.9.2-1, CMake 3.26.3 and - # * Debian 11, IntelMPI 2021, Intel Fortran 2021 - # due to malformed linker flags (the individual flags are passed as a string - # instead of a list). - # As a consequence of marking this dependency PRIVATE, every dependee of this - # target has to manually add an MPI dependency. + # MPI::MPI_Fortran cannot be made PUBLIC because the compile options of this target may cause + # errors or warnings when passed to a C compiler (e.g., when using MPICH 4.0.2 on Debian 12). + # Unfortunately we cannot manually add INTERFACE compile options because this leads to build + # problems on * JUWELS with Rocky Linux 8, ParaStationMPI 5.9.2-1, CMake 3.26.3 and * Debian 11, + # IntelMPI 2021, Intel Fortran 2021 due to malformed linker flags (the individual flags are passed + # as a string instead of a list). As a consequence of marking this dependency PRIVATE, every + # dependee of this target has to manually add an MPI dependency. target_link_libraries(negf PRIVATE MPI::MPI_Fortran) endif() @@ -129,21 +123,24 @@ if(WITH_TRANSPORT_GPU) target_link_libraries(negf PRIVATE CUDA::cudart CUDA::cusolver CUDA::cublas) endif() -if (NOT "${BLAS_LIBRARY}" STREQUAL "NONE") +if(NOT "${BLAS_LIBRARY}" STREQUAL "NONE") target_link_libraries(negf PRIVATE ${BLAS_LIBRARY}) endif() -if (NOT "${LAPACK_LIBRARY}" STREQUAL "NONE") +if(NOT "${LAPACK_LIBRARY}" STREQUAL "NONE") target_link_libraries(negf PRIVATE ${LAPACK_LIBRARY}) endif() -target_include_directories(negf PUBLIC - $ - $) +target_include_directories( + negf PUBLIC $ + $ +) -install(TARGETS negf +install( + TARGETS negf EXPORT negf-targets ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} - LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} +) if(INSTALL_INCLUDE_FILES) install(DIRECTORY ${BUILD_MOD_DIR}/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${INSTALL_MOD_DIR}) diff --git a/src/api/CMakeLists.txt b/src/api/CMakeLists.txt index 667f2222..53c9762d 100644 --- a/src/api/CMakeLists.txt +++ b/src/api/CMakeLists.txt @@ -4,43 +4,55 @@ set(all-header-deps libnegf_api.f90 bind_fortran) add_custom_target(negf-c-binding ALL DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h) set(c-header-binding-deps binding/map_negf_c.txt binding/begin_c.txt binding/end_c.txt) add_custom_command( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${all-header-deps} ${CMAKE_CURRENT_SOURCE_DIR}/${c-header-binding-deps} - COMMENT "Generating C header" - COMMAND ./bind_fortran -f c -m binding/map_negf_c.txt -b binding/begin_c.txt -n -t -e binding/end_c.txt libnegf_api.f90 > ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h + DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${all-header-deps} + ${CMAKE_CURRENT_SOURCE_DIR}/${c-header-binding-deps} + COMMENT "Generating C header" + COMMAND ./bind_fortran -f c -m binding/map_negf_c.txt -b binding/begin_c.txt -n -t -e + binding/end_c.txt libnegf_api.f90 > ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} ) list(APPEND negf-interface-headers ${CMAKE_CURRENT_BINARY_DIR}/libnegf.h) add_custom_target(negf-cpp-binding ALL DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp) -set(cpp-header-binding-deps binding/map_negf_cpp_wrapped.txt binding/begin_cpp_wrapped.txt binding/end_cpp_wrapped.txt) +set(cpp-header-binding-deps binding/map_negf_cpp_wrapped.txt binding/begin_cpp_wrapped.txt + binding/end_cpp_wrapped.txt +) add_custom_command( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${all-header-deps} ${CMAKE_CURRENT_SOURCE_DIR}/${cpp-header-binding-deps} - COMMENT "Generating C++ header" - COMMAND ./bind_fortran -f c++ -m binding/map_negf_cpp_wrapped.txt -b binding/begin_cpp_wrapped.txt -t -e binding/end_cpp_wrapped.txt libnegf_api.f90 > ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp + DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/${all-header-deps} + ${CMAKE_CURRENT_SOURCE_DIR}/${cpp-header-binding-deps} + COMMENT "Generating C++ header" + COMMAND + ./bind_fortran -f c++ -m binding/map_negf_cpp_wrapped.txt -b binding/begin_cpp_wrapped.txt -t -e + binding/end_cpp_wrapped.txt libnegf_api.f90 > ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} ) list(APPEND negf-interface-headers ${CMAKE_CURRENT_BINARY_DIR}/libnegf.hpp) -add_custom_target(negf-all-headers ALL DEPENDS negf-c-binding negf-cpp-binding ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h) +add_custom_target( + negf-all-headers ALL DEPENDS negf-c-binding negf-cpp-binding + ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h +) add_custom_command( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/lnParams.h - COMMENT "Copying lnParams.h" - COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/lnParams.h ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h + DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/lnParams.h + COMMENT "Copying lnParams.h" + COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/lnParams.h ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h ) list(APPEND negf-interface-headers ${CMAKE_CURRENT_BINARY_DIR}/lnParams.h) -add_custom_target(negf-fort-headers ALL DEPENDS negf-cpp-binding ${CMAKE_CURRENT_BINARY_DIR}/fortran.h) +add_custom_target( + negf-fort-headers ALL DEPENDS negf-cpp-binding ${CMAKE_CURRENT_BINARY_DIR}/fortran.h +) add_custom_command( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/fortran.h - DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/fortran.h - COMMENT "Copying fortran.h" - COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/fortran.h ${CMAKE_CURRENT_BINARY_DIR}/fortran.h + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/fortran.h + DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/fortran.h + COMMENT "Copying fortran.h" + COMMAND cp ${CMAKE_CURRENT_SOURCE_DIR}/fortran.h ${CMAKE_CURRENT_BINARY_DIR}/fortran.h ) list(APPEND negf-interface-headers ${CMAKE_CURRENT_BINARY_DIR}/fortran.h) if(INSTALL_INCLUDE_FILES) - install(FILES ${negf-interface-headers} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) + install(FILES ${negf-interface-headers} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) endif() diff --git a/src/contselfenergy.F90 b/src/contselfenergy.F90 index 0b3a4789..95bca63d 100644 --- a/src/contselfenergy.F90 +++ b/src/contselfenergy.F90 @@ -57,9 +57,6 @@ module ContSelfEnergy public :: SelfEnergies public :: compute_contacts public :: decimation -#:if defined("GPU") - public :: decimation_gpu -#:endif interface SelfEnergy module procedure SelfEnergy_csr @@ -76,38 +73,6 @@ module ContSelfEnergy module procedure compute_contacts_dns end interface -#:if defined("GPU") - interface - #:for PREC in PRECISIONS - #:set CTYPE = CHAR_ABBREVS['complex'][PREC] - #:set CBIND_CMPLX = ISO_C_BIND_TYPES['complex'][PREC] - #:set CBIND_REAL = ISO_C_BIND_TYPES['real'][PREC] - integer(c_int) function cu_${CTYPE}$decimation(hcublas, hcusolver, h_Go_out, h_Ao_in, h_Bo_in, h_Co_in, & - & n, tf32, ncyc, SGFACC) bind(C, name='cu_${CTYPE}$decimation') - use iso_c_binding - import cublasHandle - import cusolverDnHandle - - type(cublasHandle), value :: hcublas - type(cusolverDnHandle), value :: hcusolver - type(c_ptr), value :: h_Go_out - type(c_ptr), value :: h_Ao_in - type(c_ptr), value :: h_Bo_in - type(c_ptr), value :: h_Co_in - integer(c_size_t), value :: n - type(c_ptr), value :: ncyc - integer(c_int), value :: tf32 - real(${CBIND_REAL}$) :: SGFACC - end function - #:endfor - end interface - - interface decimation_gpu - #:for PREC in PRECISIONS - module procedure decimation_gpu_${PREC_ABBREVS[PREC]}$ - #:endfor - end interface -#:endif contains !-------------------------------------------------------------------- ! SURFACE GREEN's FUNCTION USING THE DECIMATION ITERATION @@ -209,7 +174,7 @@ subroutine surface_green(E,HC,SC,pnegf,ncyc,GS) Co=conjg(E)*SC%val(n1:n2,n3:n4)-HC%val(n1:n2,n3:n4) Co=conjg(transpose(Co)) #:if defined("GPU") - call decimation_gpu(pnegf,Go,Ao,Bo,Co,npl,.false.,ncyc) + call decimation_gpu(pnegf,Go,Ao,Bo,Co,npl,ncyc,SGFACC,loginfo=.false.) #:else call decimation(Go,Ao,Bo,Co,npl,ncyc) #:endif @@ -339,57 +304,6 @@ subroutine decimation(Go,Ao,Bo,Co,n,ncyc) end subroutine decimation -!------------------------------------------------------------------------------- -!------------------------------------------------------------------------------- -#:if defined("GPU") -#:def decimation_gpu_template(KIND, CTYPE) - - subroutine decimation_gpu_${KIND}$(negf, Go_out, Ao_in, Bo_in, Co_in, n, tf32, ncyc) - implicit none - type(Tnegf), intent(in) :: negf - integer, intent(in) :: n - complex(${KIND}$), dimension(n,n), intent(out), target :: Go_out - complex(${KIND}$), dimension(n,n), intent(in), target :: Ao_in, Bo_in, Co_in - logical, intent(in) :: tf32 - integer, intent(out), target :: ncyc - - integer :: istat, tf - type(cublasHandle) :: hh - type(cusolverDnHandle) :: hhsol - - if (tf32) then - tf = 1 - else - tf = 0 - endif - hh = negf%hcublas - hhsol = negf%hcusolver - - block - integer(c_size_t) :: n_size - - @:ASSERT(n >= 0) - n_size = int(n, kind=c_size_t) - - istat = cu_${CTYPE}$decimation(hh, hhsol, c_loc(Go_out), c_loc(Ao_in), c_loc(Bo_in), & - c_loc(Co_in), n_size, tf, c_loc(ncyc), & - real(SGFACC,${KIND}$)*n*n) - end block - - end subroutine decimation_gpu_${KIND}$ - -#:enddef decimation_gpu_template - -#:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] - #:set CTYPE = CHAR_ABBREVS['complex'][PREC] - - $:decimation_gpu_template(KIND,CTYPE) -#:endfor - - -#:endif - !------------------------------------------------------------------------------- subroutine compute_contacts_csr(Ec,pnegf,avncyc,Tlc,Tcl,SelfEneR,GS) diff --git a/src/cudautils.F90 b/src/cudautils.F90 index 984b23ce..eaeb480c 100644 --- a/src/cudautils.F90 +++ b/src/cudautils.F90 @@ -28,7 +28,8 @@ module cudautils use ln_precision use mat_def use lib_param - use iso_c_binding, only : c_int, c_loc, c_null_ptr, c_ptr, c_size_t + use iso_c_binding, only : c_int, c_int32_t, c_int64_t, c_size_t, & + & c_loc, c_null_ptr, c_ptr, c_bool use, intrinsic :: ieee_arithmetic use ln_messages implicit none @@ -61,6 +62,7 @@ module cudautils public :: copy_mat_gpu public :: asum_gpu public :: dagger_gpu + public :: decimation_gpu public :: checksum @@ -72,7 +74,7 @@ module cudautils public :: getDeviceCount public :: getDevice public :: setDevice - !public :: getDeviceProperties + public :: getDeviceProperties public :: getDevMemInfo @@ -220,6 +222,12 @@ module cudautils #:endfor end interface dagger_gpu + interface decimation_gpu + #:for PREC in PRECISIONS + module procedure decimation_gpu_${PREC_ABBREVS[PREC]}$ + #:endfor + end interface decimation_gpu + interface checksum #:for PREC in PRECISIONS module procedure checksum_${PREC_ABBREVS[PREC]}$ @@ -236,28 +244,24 @@ module cudautils interface integer(c_int) function cu_cublasInit(hcublas) bind(C, name='cu_cublasInit') use iso_c_binding - import cublasHandle - type(cublasHandle) :: hcublas + type(c_ptr) :: hcublas end function cu_cublasInit integer(c_int) function cu_cublasFinalize(hcublas) bind(C, name='cu_cublasFinalize') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas + type(c_ptr), value :: hcublas end function cu_cublasFinalize integer(c_int) function cu_cusolverInit(hcusolver) bind(C, name='cu_cusolverInit') use iso_c_binding - import cusolverDnHandle - type(cusolverDnHandle) :: hcusolver + type(c_ptr) :: hcusolver end function cu_cusolverInit integer(c_int) function cu_cusolverFinalize(hcusolver) bind(C, name='cu_cusolverFinalize') use iso_c_binding - import cusolverDnHandle - type(cusolverDnHandle), value :: hcusolver + type(c_ptr), value :: hcusolver end function cu_cusolverFinalize - + integer(c_int) function cudaDeviceSynchronize() bind(C, name='cudaDeviceSynchronize') use iso_c_binding end function cudaDeviceSynchronize @@ -273,18 +277,17 @@ integer(c_int) function cu_cudaGetDeviceCount(count) bind(C, name='cu_cudaGetDev integer(c_int) :: count end function cu_cudaGetDeviceCount - integer(c_int) function cu_cudaGetDeviceProperties(device,prop) & + integer(c_int) function cu_cudaGetDeviceProperties(device) & & bind(C, name='cu_cudaGetDeviceProperties') use iso_c_binding integer(c_int) :: device - type(c_ptr) :: prop end function cu_cudaGetDeviceProperties - integer(c_int) function cu_cudaMallocAsync(d_A, siz) bind(C, name='cu_cudaMallocAsync') + integer(c_int) function cu_cudaMallocAsync(d_A, msize) bind(C, name='cu_cudaMallocAsync') use iso_c_binding ! not by value since pointer has to be initialized, hence pass its reference type(c_ptr) :: d_A - integer(c_size_t), value :: siz + integer(c_size_t), value :: msize end function cu_cudaMallocAsync integer(c_int) function cu_cudaGetDevice(device) bind(C, name='cu_cudaGetDevice') @@ -345,30 +348,26 @@ end function cu_meminfo #:set C_BIND_CMPLX = ISO_C_BIND_TYPES['complex'][PREC] #:set C_BIND_REAL = ISO_C_BIND_TYPES['real'][PREC] ! C = alpha*A*B + beta*C - integer(c_int) function cu_${CTYPE}$multMat(hcublas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger) & - & bind(C, name='cu_${CTYPE}$multMat') + integer(c_int) function cu_${CTYPE}$matmul(hcublas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger) & + & bind(C, name='cu_${CTYPE}$matmul') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas - integer(c_size_t), value :: m - integer(c_size_t), value :: n - integer(c_size_t), value :: k + type(c_ptr), value :: hcublas + integer(c_int32_t), value :: m + integer(c_int32_t), value :: n + integer(c_int32_t), value :: k complex(${C_BIND_CMPLX}$) :: alpha type(c_ptr), value :: d_A type(c_ptr), value :: d_B complex(${C_BIND_CMPLX}$) :: beta type(c_ptr), value :: d_C integer(c_int), value :: dagger - end function cu_${CTYPE}$multMat + end function cu_${CTYPE}$matmul - integer(c_int) function cu_${CTYPE}$inverse(hcublas, hcusolver, d_A, d_Ainv, N) & + integer(c_int) function cu_${CTYPE}$inverse(hcusolver, d_A, d_Ainv, N) & & bind(C, name='cu_${CTYPE}$inverse') use iso_c_binding - import cublasHandle - import cusolverDnHandle - type(cusolverDnHandle), value :: hcusolver - type(cublasHandle), value :: hcublas - integer(c_size_t), value :: N + type(c_ptr), value :: hcusolver + integer(c_int32_t), value :: N type(c_ptr), value :: d_A type(c_ptr), value :: d_Ainv end function cu_${CTYPE}$inverse @@ -376,10 +375,9 @@ end function cu_${CTYPE}$inverse integer(c_int) function cu_${CTYPE}$matsum(hcublas, m, n, alpha, d_A, beta, d_B, d_C, dagger) & & bind(C, name='cu_${CTYPE}$matsum') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas - integer(c_size_t), value :: m - integer(c_size_t), value :: n + type(c_ptr), value :: hcublas + integer(c_int32_t), value :: m + integer(c_int32_t), value :: n integer(c_int), value :: dagger complex(${C_BIND_CMPLX}$) :: alpha type(c_ptr), value :: d_A @@ -392,40 +390,56 @@ integer(c_int) function cu_${CTYPE}$initmat(d_A, nrow) & & bind(C, name='cu_${CTYPE}$initmat') use iso_c_binding type(c_ptr), value :: d_A - integer(c_size_t), value :: nrow + integer(c_int32_t), value :: nrow end function cu_${CTYPE}$initmat real(${C_BIND_REAL}$) function cu_${CTYPE}$trace(hcublas, d_A, nrow, h_tun, mask_present) & & bind(C, name='cu_${CTYPE}$trace') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas + type(c_ptr), value :: hcublas type(c_ptr), value :: d_A type(c_ptr), value :: h_tun - integer(c_size_t), value :: nrow + integer(c_int32_t), value :: nrow integer(c_int), value :: mask_present end function cu_${CTYPE}$trace integer(c_int) function cu_${CTYPE}$matcopy(hcublas, d_A, d_B, msize) & & bind(C, name='cu_${CTYPE}$matcopy') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas + type(c_ptr), value :: hcublas type(c_ptr), value :: d_A type(c_ptr), value :: d_B integer(c_size_t), value :: msize end function cu_${CTYPE}$matcopy - integer(c_int) function cu_${CTYPE}$asum(hcublas, d_A, summ, N) & + integer(c_int) function cu_${CTYPE}$asum(hcublas, d_A, summ, nelems) & & bind(C, name='cu_${CTYPE}$asum') use iso_c_binding - import cublasHandle - type(cublasHandle), value :: hcublas + type(c_ptr), value :: hcublas type(c_ptr), value :: d_A real(${C_BIND_REAL}$) :: summ - integer(c_size_t), value :: N + integer(c_size_t), value :: nelems end function cu_${CTYPE}$asum + integer(c_int) function cu_${CTYPE}$decimation(hcublas, hcusolver, h_Go_out, & + & h_Ao_in, h_Bo_in, h_Co_in, & + & n, ncyc, accuracy, loginfo) bind(C, name='cu_${CTYPE}$decimation') + use iso_c_binding + import cublasHandle + import cusolverDnHandle + + type(cublasHandle), value :: hcublas + type(cusolverDnHandle), value :: hcusolver + type(c_ptr), value :: h_Go_out + type(c_ptr), value :: h_Ao_in + type(c_ptr), value :: h_Bo_in + type(c_ptr), value :: h_Co_in + integer(c_int32_t), value :: n + type(c_ptr), value :: ncyc + real(${C_BIND_REAL}$), value :: accuracy + logical(kind=c_bool), value :: loginfo + end function + #:endfor end interface @@ -437,25 +451,25 @@ end function cu_${CTYPE}$asum subroutine cublasInitialize(hcublas) type(cublasHandle), intent(inout) :: hcublas integer :: err - err = cu_cublasInit(hcublas) + err = cu_cublasInit(hcublas%ptr) end subroutine cublasInitialize subroutine cublasFinalize(hcublas) type(cublasHandle), intent(inout) :: hcublas integer :: err - err = cu_cublasFinalize(hcublas) + err = cu_cublasFinalize(hcublas%ptr) end subroutine cublasFinalize subroutine cusolverInitialize(hcusolver) type(cusolverDnHandle), intent(inout) :: hcusolver integer :: err - err = cu_cusolverInit(hcusolver) + err = cu_cusolverInit(hcusolver%ptr) end subroutine cusolverInitialize subroutine cusolverFinalize(hcusolver) type(cusolverDnHandle), intent(inout) :: hcusolver integer :: err - err = cu_cusolverFinalize(hcusolver) + err = cu_cusolverFinalize(hcusolver%ptr) end subroutine cusolverFinalize subroutine getDeviceCount(count) @@ -482,6 +496,12 @@ subroutine getDevMemInfo(freemem, totalmem) err = cu_meminfo(freemem, totalmem) end subroutine getDevMemInfo + subroutine getDeviceProperties(dev) + integer(4) :: dev + integer :: err + err = cu_cudaGetDeviceProperties(dev) + end subroutine getDeviceProperties + subroutine waitForGPU() integer :: err err = cudaDeviceSynchronize() @@ -495,7 +515,9 @@ subroutine createGPU_${KIND}$(A) type(${MTYPE}$), intent(in) :: A integer :: err - err = cu_createMat(A%d_addr, size(A%val,kind=c_size_t)*${CUDATYPE}$) + integer(c_size_t) :: bytecount + bytecount = size(A%val,kind=c_size_t)*${CUDATYPE}$ + err = cu_createMat(A%d_addr, bytecount) end subroutine createGPU_${KIND}$ #:enddef createGPU_template @@ -504,7 +526,9 @@ subroutine createGPU_async_${KIND}$(A) use iso_c_binding, only: c_size_t type(${MTYPE}$), intent(in) :: A integer :: err - err = cu_cudaMallocAsync(A%d_addr, size(A%val,kind=c_size_t)*${CUDATYPE}$) + integer(c_size_t) :: bytecount + bytecount = size(A%val,kind=c_size_t)*${CUDATYPE}$ + err = cu_cudaMallocAsync(A%d_addr, bytecount) end subroutine createGPU_async_${KIND}$ #:enddef createGPU_async_template @@ -515,9 +539,11 @@ subroutine createGPU_only_async_${KIND}$(A, m, n) integer, intent(in) :: m, n integer :: err + integer(c_size_t) :: bytecount A%nrow = m A%ncol = n - err = cu_cudaMallocAsync(A%d_addr, int(m, kind=c_size_t) * n *${CUDATYPE}$) + bytecount = int(m,kind=c_size_t)*int(n,kind=c_size_t)*${CUDATYPE}$ + err = cu_cudaMallocAsync(A%d_addr, bytecount) @:ASSERT(err == 0) end subroutine createGPU_only_async_${KIND}$ #:enddef createGPU_only_async_template @@ -527,8 +553,13 @@ end subroutine createGPU_only_async_${KIND}$ subroutine copyToGPU_${KIND}$(A) type(${MTYPE}$), intent(in), target :: A integer :: err - !call createGPU(A) - err = cu_copyMatH2D(c_loc(A%val), A%d_addr, size(A%val,kind=c_size_t)*${CUDATYPE}$) + integer(c_size_t) :: bytecount + bytecount = size(A%val,kind=c_size_t)*${CUDATYPE}$ + err = cu_copyMatH2D(c_loc(A%val(1,1)), A%d_addr, bytecount) + if (err /= 0) then + print*, 'COPY ',c_loc(A%val(1,1)), 'TO ', A%d_addr, bytecount + error stop 'COPY ERROR !!' + end if end subroutine copyToGPU_${KIND}$ #:enddef copyToGPU_template @@ -538,7 +569,13 @@ subroutine copyToGPU_async_${KIND}$(A) use iso_c_binding, only: c_size_t type(${MTYPE}$), intent(in), target :: A integer :: err - err = cu_copyMatH2D_async(c_loc(A%val), A%d_addr, size(A%val,kind=c_size_t)*${CUDATYPE}$) + integer(c_size_t) :: bytecount + bytecount = size(A%val,kind=c_size_t)*${CUDATYPE}$ + err = cu_copyMatH2D_async(c_loc(A%val(1,1)), A%d_addr, bytecount) + if (err /= 0) then + print*, 'COPY ',c_loc(A%val(1,1)), 'TO ', A%d_addr, bytecount + error stop 'COPY ERROR !!' + end if end subroutine copyToGPU_async_${KIND}$ #:enddef copyToGPU_async_template @@ -547,7 +584,9 @@ end subroutine copyToGPU_async_${KIND}$ subroutine copyFromGPU_${KIND}$(A) type(${MTYPE}$), intent(in), target :: A integer :: err - err = cu_copyMatD2H(c_loc(A%val), A%d_addr, size(A%val,kind=c_size_t)*${CUDATYPE}$) + integer(c_size_t) :: bytecount + bytecount = size(A%val,kind=c_size_t)*${CUDATYPE}$ + err = cu_copyMatD2H(c_loc(A%val(1,1)), A%d_addr, bytecount) end subroutine copyFromGPU_${KIND}$ #:enddef copyFromGPU_template @@ -602,31 +641,31 @@ subroutine matmul_gpu_${KIND}$(hcublas, alpha, A, B, beta, C, dagger) character(*), intent(in), optional :: dagger integer :: istat - integer(c_size_t) :: ncol_A - integer(c_size_t) :: ncol_C - integer(c_size_t) :: nrow_B - integer(c_size_t) :: nrow_C + integer(c_int32_t) :: ncol_A + integer(c_int32_t) :: ncol_C + integer(c_int32_t) :: nrow_B + integer(c_int32_t) :: nrow_C @:ASSERT(A%ncol >= 0) - ncol_A = int(A%ncol, kind=c_size_t) + ncol_A = int(A%ncol, kind=c_int32_t) @:ASSERT(C%ncol >= 0) - ncol_C = int(C%ncol, kind=c_size_t) + ncol_C = int(C%ncol, kind=c_int32_t) @:ASSERT(B%nrow >= 0) - nrow_B = int(B%nrow, kind=c_size_t) + nrow_B = int(B%nrow, kind=c_int32_t) @:ASSERT(C%nrow >= 0) - nrow_C = int(C%nrow, kind=c_size_t) + nrow_C = int(C%nrow, kind=c_int32_t) if (.not.present(dagger)) then - istat = cu_${CTYPE}$multMat(hcublas, nrow_C, ncol_C, ncol_A, alpha, A%d_addr, & + istat = cu_${CTYPE}$matmul(hcublas%ptr, nrow_C, ncol_C, ncol_A, alpha, A%d_addr, & & B%d_addr, beta, C%d_addr, 0) else select case(dagger) case('dag_1st') - istat = cu_${CTYPE}$multMat(hcublas, nrow_C, ncol_C, nrow_B, alpha, A%d_addr, & + istat = cu_${CTYPE}$matmul(hcublas%ptr, nrow_C, ncol_C, nrow_B, alpha, A%d_addr, & & B%d_addr, beta, C%d_addr, 1) case('dag_2nd') - istat = cu_${CTYPE}$multMat(hcublas, nrow_C, ncol_C, ncol_A, alpha, A%d_addr, & + istat = cu_${CTYPE}$matmul(hcublas%ptr, nrow_C, ncol_C, ncol_A, alpha, A%d_addr, & & B%d_addr, beta, C%d_addr, 2) case default call error_msg('Error in matmul_gpu') @@ -645,13 +684,13 @@ subroutine inverse_gpu_${KIND}$(hcublas, hcusolver, A, Ainv, err) integer, intent(out) :: err integer :: istat - integer(c_size_t) :: nrow + integer(c_int32_t) :: nrow call init_gpu_${KIND}$(Ainv) @:ASSERT(A%nrow >= 0) - nrow = int(A%nrow, kind=c_size_t) - istat = cu_${CTYPE}$inverse(hcublas,hcusolver, A%d_addr, Ainv%d_addr, nrow) + nrow = int(A%nrow, kind=c_int32_t) + istat = cu_${CTYPE}$inverse(hcusolver%ptr, A%d_addr, Ainv%d_addr, nrow) err = istat end subroutine inverse_gpu_${KIND}$ @@ -666,9 +705,10 @@ subroutine kernelsum_gpu_${KIND}$(C,alpha,A,beta,B) complex(${KIND}$), intent(in) :: beta integer :: istat + integer(c_size_t) :: nelems - istat = cu_${CTYPE}$kernelsum(C%d_addr, alpha, A%d_addr, beta, B%d_addr, & - & size(A%val, kind=c_size_t)) + nelems = int(A%val,kind=c_size_t) + istat = cu_${CTYPE}$kernelsum(C%d_addr, alpha, A%d_addr, beta, B%d_addr, nelems) end subroutine kernelsum_gpu_${KIND}$ #:enddef kernelsum_gpu_template @@ -684,24 +724,24 @@ subroutine matsum_gpu_${KIND}$(hcublas, alpha, A, beta, B, C, dagger) character(*), intent(in), optional :: dagger integer :: istat - integer(c_size_t) :: ncol_C - integer(c_size_t) :: nrow_C + integer(c_int32_t) :: ncol_C + integer(c_int32_t) :: nrow_C @:ASSERT(C%ncol >= 0) - ncol_C = int(C%ncol, kind=c_size_t) + ncol_C = int(C%ncol, kind=c_int32_t) @:ASSERT(C%nrow >= 0) - nrow_C = int(C%nrow, kind=c_size_t) + nrow_C = int(C%nrow, kind=c_int32_t) if (.not.present(dagger)) then - istat = cu_${CTYPE}$matsum(hcublas, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & + istat = cu_${CTYPE}$matsum(hcublas%ptr, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & & C%d_addr, 0) else if (dagger == 'dag_1st') then - istat = cu_${CTYPE}$matsum(hcublas, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & + istat = cu_${CTYPE}$matsum(hcublas%ptr, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & & C%d_addr, 1) endif if (dagger == 'dag_2nd') then - istat = cu_${CTYPE}$matsum(hcublas, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & + istat = cu_${CTYPE}$matsum(hcublas%ptr, nrow_C, ncol_C, alpha, A%d_addr, beta, B%d_addr, & & C%d_addr, 2) endif endif @@ -730,10 +770,10 @@ subroutine init_gpu_${KIND}$(A) type(${MTYPE}$), intent(inout) :: A integer :: istat - integer(c_size_t) :: nrow + integer(c_int32_t) :: nrow @:ASSERT(A%nrow >= 0) - nrow = int(A%nrow, kind=c_size_t) + nrow = int(A%nrow, kind=c_int32_t) istat = cu_${CTYPE}$initmat(A%d_addr, nrow) end subroutine init_gpu_${KIND}$ @@ -747,16 +787,16 @@ subroutine trace_gpu_${KIND}$(hcublas, A, trace, tun_mask) logical, intent(in), optional, target :: tun_mask(:) type(c_ptr) :: dummy - integer(c_size_t) :: nrow + integer(c_int32_t) :: nrow dummy = c_null_ptr @:ASSERT(A%nrow >= 0) - nrow = int(A%nrow, kind=c_size_t) + nrow = int(A%nrow, kind=c_int32_t) if (.not.present(tun_mask)) then - trace = cu_${CTYPE}$trace(hcublas, A%d_addr, nrow, dummy, 0) + trace = cu_${CTYPE}$trace(hcublas%ptr, A%d_addr, nrow, dummy, 0) else - trace = cu_${CTYPE}$trace(hcublas, A%d_addr, nrow, c_loc(tun_mask), 1) + trace = cu_${CTYPE}$trace(hcublas%ptr, A%d_addr, nrow, c_loc(tun_mask), 1) endif end subroutine trace_gpu_${KIND}$ @@ -769,12 +809,12 @@ subroutine copy_mat_gpu_${KIND}$(hcublas, A, Acopy) type(${MTYPE}$), intent(inout) :: Acopy integer :: istat - integer(c_size_t) :: n + integer(c_size_t) :: nelems @:ASSERT(A%ncol >= 0) @:ASSERT(A%nrow >= 0) - n = int(A%ncol, kind=c_size_t) * int(A%nrow, kind=c_size_t) - istat = cu_${CTYPE}$matcopy(hcublas, A%d_addr, Acopy%d_addr, n) + nelems = int(A%ncol, kind=c_size_t) * int(A%nrow, kind=c_size_t) + istat = cu_${CTYPE}$matcopy(hcublas%ptr, A%d_addr, Acopy%d_addr, nelems) end subroutine copy_mat_gpu_${KIND}$ #:enddef copy_mat_gpu_template @@ -786,12 +826,12 @@ subroutine asum_gpu_${KIND}$(hcublas, A, summ) real(${KIND}$), intent(out) :: summ integer :: istat - integer(c_size_t) :: n + integer(c_size_t) :: nelems @:ASSERT(A%ncol >= 0) @:ASSERT(A%nrow >= 0) - n = int(A%ncol, kind=c_size_t) * int(A%nrow, kind=c_size_t) - istat = cu_${CTYPE}$asum(hcublas, A%d_addr, summ, n) + nelems = int(A%ncol, kind=c_size_t) * int(A%nrow, kind=c_size_t) + istat = cu_${CTYPE}$asum(hcublas%ptr, A%d_addr, summ, nelems) end subroutine asum_gpu_${KIND}$ #:enddef asum_gpu_template @@ -909,6 +949,39 @@ subroutine checksum_${KIND}$(hcublas, A, nome) end subroutine checksum_${KIND}$ #:enddef checksum_template +#:def decimation_gpu_template(KIND, CTYPE) + subroutine decimation_gpu_${KIND}$(negf, Go_out, Ao_in, Bo_in, Co_in, n, ncyc, accuracy, loginfo) + implicit none + type(Tnegf), intent(in) :: negf + integer, intent(in) :: n + complex(${KIND}$), dimension(n,n), intent(out), target :: Go_out + complex(${KIND}$), dimension(n,n), intent(in), target :: Ao_in, Bo_in, Co_in + integer, intent(out), target :: ncyc + real(${KIND}$), intent(in) :: accuracy + logical, intent(in) :: loginfo + + integer :: istat + type(cublasHandle) :: hh + type(cusolverDnHandle) :: hhsol + + hh = negf%hcublas + hhsol = negf%hcusolver + + block + integer(c_int32_t) :: nrows + logical(kind=c_bool) :: logg + + @:ASSERT(n >= 0) + nrows = int(n, kind=c_int32_t) + logg = logical(loginfo, kind=c_bool) + + istat = cu_${CTYPE}$decimation(hh, hhsol, c_loc(Go_out), c_loc(Ao_in), c_loc(Bo_in), & + c_loc(Co_in), nrows, c_loc(ncyc), accuracy, logg) + end block + + end subroutine decimation_gpu_${KIND}$ + +#:enddef decimation_gpu_template #:for PREC in PRECISIONS #:set KIND = PREC_ABBREVS[PREC] @@ -966,5 +1039,6 @@ end subroutine checksum_${KIND}$ $:checksum_template(KIND, CTYPE, MTYPE, CUDATYPE) + $:decimation_gpu_template(KIND, CTYPE) #:endfor end module cudautils diff --git a/src/elphinel.F90 b/src/elphinel.F90 index d2abc65d..dd4dce38 100644 --- a/src/elphinel.F90 +++ b/src/elphinel.F90 @@ -234,7 +234,7 @@ subroutine ElPhonNonPO_init(this, comm, struct, basis, coupling, wq, Temp, & logical, intent(in) :: tridiag this%descriptor = & - & "Electron-Phonon inelastic model for polar-optical phonons" + & "Electron-Phonon inelastic model for non-PO phonons" this%cart_comm = comm this%scba_niter = niter diff --git a/src/gpuroutines.cu b/src/gpuroutines.cu index 6b6b1777..dadda7d8 100644 --- a/src/gpuroutines.cu +++ b/src/gpuroutines.cu @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -58,6 +59,19 @@ } \ } while(false) #endif +#define CUDA_CHECK(error) \ + do { \ + if(!(error == cudaSuccess)) { \ + std::fprintf( \ + stderr, "CUDA error in %s line %u \n", __FILE__, \ + __LINE__ \ + ); \ + std::fprintf( \ + stderr, "ERROR \n %s \n", cudaGetErrorString(error) \ + ); \ + std::exit(EXIT_FAILURE); \ + } \ + } while(false) #define LIBNEGF_CHECK_INFO(info_var) \ @@ -118,15 +132,15 @@ struct get_real { template __global__ void addKernel( Number* c, Number alpha, const Number* a, Number beta, const Number* b, - size_t size + size_t num_elements ) { ENFORCE(c); ENFORCE(a); ENFORCE(b); - auto i = blockIdx.x * blockDim.x + threadIdx.x; + size_t i = blockIdx.x * blockDim.x + threadIdx.x; - if(i < size) { + if(i < num_elements) { c[i].x = (alpha.x * a[i].x - alpha.y * a[i].y) + (beta.x * b[i].x - beta.y * b[i].y); c[i].y = (alpha.x * a[i].y + alpha.y * a[i].x) + @@ -165,17 +179,18 @@ __global__ void hermitian(cuComplex *odata, const cuComplex *idata) * Initializes a square complex matrix as the identity matrix. */ template -__global__ void initKernel(Number* a, size_t nrow) { +__global__ void initKernel(Number* a, int32_t nrow) { using Real = typename get_real::type; ENFORCE(a); - auto size = nrow * nrow; - auto i = blockDim.x * blockIdx.x + threadIdx.x; + size_t n64 = static_cast(nrow); + size_t num_elements = n64 * n64; + size_t i = blockDim.x * blockIdx.x + threadIdx.x; auto one = Number{Real{1}}; auto zero = Number{Real{0}}; - if(i < size) { + if(i < num_elements) { if(i % (nrow + 1) == 0) { a[i] = one; } else { @@ -186,7 +201,7 @@ __global__ void initKernel(Number* a, size_t nrow) { template -__global__ void initArrayWithOnes(Number* a, size_t nrow) { +__global__ void initArrayWithOnes(Number* a, int32_t nrow) { static_assert(std::is_same::type>::value); assert(a); @@ -198,49 +213,119 @@ __global__ void initArrayWithOnes(Number* a, size_t nrow) { } } - -/** - * Computes the trace of a matrix A. - * - * Optionally, a bit mask can be passed to the function. If a mask is present, - * then a diagonal element (i,i) is only considered for the trace computation - * if the i-th value in mask is nonzero. - * - * @param[in] mask A bitmask of length nrow indicating which rows to ignore - * (zero means ignore). +/* *************************************************************** + * Extract diagonal from a square matrix + * *************************************************************** */ + template::type> -__global__ void -traceKernel(Number* a, size_t nrow, Real* trace, bool* mask, int mask_present) { - ENFORCE(a); - ENFORCE(trace); +__global__ void getdiagKernel(const Number* A, int n, Real* out, + const bool* mask, int mask_present) +{ + ENFORCE(A); + ENFORCE(out); + + int i = blockIdx.x * blockDim.x + threadIdx.x; + + if (i < n) { + if (mask_present == 0 || mask[i]) { + out[i] = A[i * n + i].x; + } else { + out[i] = 0; + } + } +} + +/* *************************************************************** + * Simple trace kernel with a mask + * *************************************************************** + */ +/* +template +__global__ void traceKernel(const Number* A, int32_t nrow, Number* result, + bool* mask, int mask_present) +{ + ENFORCE(A); + ENFORCE(result); ENFORCE(mask || mask_present == 0); ENFORCE(mask_present == 0 || mask_present == 1); - auto size = nrow * nrow; - auto i = blockDim.x * blockIdx.x + threadIdx.x; + size_t i = blockIdx.x * blockDim.x + threadIdx.x; if(mask_present == 0) { - if(i < size) { - if(i % (nrow + 1) == 0) { - trace[i % nrow] = a[i].x; - } - } + if (i < nrow) { + Number val = A[i * nrow + i]; + atomicAdd(&((*result).x), val.x); + atomicAdd(&((*result).y), val.y); + } } - if(mask_present == 1) { - if(i < size) { - if(i % (nrow + 1) == 0) { - if(mask[i % nrow]) { - trace[i % nrow] = a[i].x; - } else { - trace[i % nrow] = 0.0; - } - } - } + else + { + if (i < nrow && mask[i]) { + Number val = A[i * nrow + i]; + atomicAdd(&((*result).x), val.x); + atomicAdd(&((*result).y), val.y); + } } } +*/ + +__device__ inline float absval(cuComplex z) { return cuCabsf(z); } +__device__ inline double absval(cuDoubleComplex z){ return cuCabs(z); } + +__device__ __forceinline__ float atomicMaxAbsFloat(float* addr, float value) { + // assumes value >= 0 + return __int_as_float(atomicMax((int*)addr, __float_as_int(value))); +} + +__device__ __forceinline__ double atomicMaxAbsDouble(double* addr, double value) +{ + // assumes value >= 0 + return __longlong_as_double( + atomicMax((unsigned long long*)addr, __double_as_longlong(value)) + ); +} +/* ------------------------------------------------------------------------ + * Kernel to find maxval(abs(A)) + * ------------------------------------------------------------------------ + */ +template::type> +__global__ void maxabs_kernel(const Number* __restrict__ data, + size_t n, + Real* __restrict__ norm) +{ + extern __shared__ unsigned char smem[]; + Real* sdata = reinterpret_cast(smem); + + unsigned int tid = threadIdx.x; + size_t i = blockDim.x * blockIdx.x + tid; + + Real val = 0; + if (i < n) val = absval(data[i]); + + sdata[tid] = val; + + __syncthreads(); + + // reduction to find block maximum + for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) + { + if (tid < s) + sdata[tid] = fmax(sdata[tid], sdata[tid + s]); + __syncthreads(); + } + + // final block max reduced atomically to global result + if (tid == 0) { + if constexpr (std::is_same::value) + atomicMaxAbsFloat(norm, sdata[0]); + else + atomicMaxAbsDouble(norm, sdata[0]); + } +} + /* * DATA MOVEMENT ROUTINES */ @@ -248,7 +333,7 @@ traceKernel(Number* a, size_t nrow, Real* trace, bool* mask, int mask_present) { extern "C" int cu_createMat(void** d_A, size_t bytecount) { ENFORCE(d_A); cudaError_t err = cudaMalloc(d_A, bytecount); - //printf("create mat at GPU Address: %p \n",*d_A); + CUDA_CHECK(err); return err; } @@ -264,23 +349,21 @@ extern "C" int cu_cudaFreeAsync(void** d_A) { extern "C" int cu_cudaMallocAsync(void** d_A, size_t bytecount) { ENFORCE(d_A); cudaError_t err = cudaMallocAsync(d_A, bytecount, 0); - //printf("create mat at GPU Address: %p \n",*d_A); - ENFORCE(err == cudaSuccess); + CUDA_CHECK(err); return err; } extern "C" int cu_copyMatH2D(void* h_A, void* d_A, size_t bytecount) { ENFORCE(h_A); ENFORCE(d_A); - // printf("copy %p to %p\n",h_A,d_A); cudaError_t err = cudaMemcpy(d_A, h_A, bytecount, cudaMemcpyHostToDevice); + CUDA_CHECK(err); return err; } extern "C" int cu_copyMatH2D_async(void* h_A, void* d_A, size_t bytecount) { ENFORCE(h_A); ENFORCE(d_A); - // printf("copy %p to %p\n",h_A,d_A); cudaError_t err = cudaMemcpyAsync(d_A, h_A, bytecount, cudaMemcpyHostToDevice); return err; @@ -289,8 +372,24 @@ extern "C" int cu_copyMatH2D_async(void* h_A, void* d_A, size_t bytecount) { extern "C" int cu_copyMatD2H(void* h_A, void* d_A, size_t bytecount) { ENFORCE(h_A); ENFORCE(d_A); - cudaError_t err = cudaMemcpy(h_A, d_A, bytecount, cudaMemcpyDeviceToHost); + CUDA_CHECK(err); + return err; +} + +extern "C" int cu_copyMatD2D_async(void* d_A, void* d_A_in, size_t bytecount) { + ENFORCE(d_A); + ENFORCE(d_A_in); + cudaError_t err = + cudaMemcpyAsync(d_A, d_A_in, bytecount, cudaMemcpyDeviceToDevice); + return err; +} + +extern "C" int cu_copyMatD2D(void* d_A, void* d_A_in, size_t bytecount) { + ENFORCE(d_A); + ENFORCE(d_A_in); + cudaError_t err = cudaMemcpy(d_A, d_A_in, bytecount, cudaMemcpyDeviceToDevice); + CUDA_CHECK(err); return err; } @@ -311,14 +410,14 @@ extern "C" int cu_deleteMat(void** d_A) { extern "C" int cu_cudaGetDeviceCount(int* count) { ENFORCE(count); cudaError_t err = cudaGetDeviceCount(count); - ENFORCE(err == cudaSuccess); + CUDA_CHECK(err); return err; } extern "C" int cu_cudaGetDeviceProperties(int device) { cudaDeviceProp prop; cudaError_t err = cudaGetDeviceProperties(&prop, device); - ENFORCE(err == cudaSuccess); + CUDA_CHECK(err); printf(" Found GPU: Device Name: %s\n", prop.name); printf(" TotalMemory: %lu\n", (unsigned long)prop.totalGlobalMem); @@ -329,45 +428,54 @@ extern "C" int cu_cudaGetDeviceProperties(int device) { extern "C" int cu_cudaGetDevice(int* device) { cudaError_t err = cudaGetDevice(device); - assert(err == cudaSuccess); + CUDA_CHECK(err); return err; } extern "C" int cu_cudaSetDevice(int count) { cudaError_t err = cudaSetDevice(count); - ENFORCE(err == cudaSuccess); + CUDA_CHECK(err); return err; } -extern "C" int cu_cublasInit(cublasHandle_t* hcublas) { - ENFORCE(hcublas); - cublasStatus_t err = cublasCreate(hcublas); +extern "C" int cu_cublasInit(void** hpblas) { + ENFORCE(hpblas); + + cublasHandle_t hb; + cublasStatus_t err = cublasCreate(&hb); ENFORCE(err == CUBLAS_STATUS_SUCCESS); if(err != CUBLAS_STATUS_SUCCESS) { printf("cublas create error: %d\n", err); } // printf("hcublas Addr: %p \n",*hcublas); + *hpblas = hb; return err; } -extern "C" int cu_cublasFinalize(cublasHandle_t hcublas) { +extern "C" int cu_cublasFinalize(void *hpblas) { + ENFORCE(hpblas); + cublasHandle_t hcublas = static_cast(hpblas); cublasStatus_t err = cublasDestroy(hcublas); ENFORCE(err == CUBLAS_STATUS_SUCCESS); return err; } -extern "C" int cu_cusolverInit(cusolverDnHandle_t* hcusolver) { - ENFORCE(hcusolver); - cusolverStatus_t err = cusolverDnCreate(hcusolver); +extern "C" int cu_cusolverInit(void** hpsolver) { + ENFORCE(hpsolver); + cusolverDnHandle_t hs; + cusolverStatus_t err = cusolverDnCreate(&hs); ENFORCE(err == CUSOLVER_STATUS_SUCCESS); if(err != 0) { printf("cusolver create error: %d\n", err); } // printf("hcusolver Addr: %p \n",*hcusolver); + *hpsolver = hs; return err; } -extern "C" int cu_cusolverFinalize(cusolverDnHandle_t hcusolver) { +extern "C" int cu_cusolverFinalize(void *hpsolver) { + ENFORCE(hpsolver); + cusolverDnHandle_t hcusolver = static_cast(hpsolver); cusolverStatus_t err = cusolverDnDestroy(hcusolver); ENFORCE(err == CUSOLVER_STATUS_SUCCESS); return err; @@ -383,9 +491,9 @@ extern "C" int cu_cusolverFinalize(cusolverDnHandle_t hcusolver) { * * op(.) indicates if the matrix or its complex-conjugate is used. The allowed * values for dagger are: - * * `dagger == 0`: compute C ≔ α A · B + β C - * * `dagger == 1`: compute C ≔ α A^* · B + β C - * * `dagger == 2`: compute C ≔ α A · B^* + β C + * * `dagger == 0`: compute C = α A · B + β C + * * `dagger == 1`: compute C = α A^* · B + β C + * * `dagger == 2`: compute C = α A · B^* + β C * * @param[in] m The number of rows of C and op(A). * @param[in] n The number of columns of C and op(B). @@ -393,8 +501,8 @@ extern "C" int cu_cusolverFinalize(cusolverDnHandle_t hcusolver) { * @param[in] dagger A shorthand for various combinations of op(A), op(B). */ template -int cu_multMat( - cublasHandle_t hcublas, size_t m, size_t n, size_t k, const Number* alpha, +int cu_matmul( + void *hpblas, int32_t m, int32_t n, int32_t k, const Number* alpha, const Number* d_A, const Number* d_B, const Number* beta, Number* d_C, int dagger ) { @@ -404,6 +512,7 @@ int cu_multMat( ENFORCE(beta); ENFORCE(dagger == 0 || dagger == 1 || dagger == 2); + cublasHandle_t hcublas = static_cast(hpblas); cublasStatus_t err; if(dagger == 0) { err = libnegf::cublasGemm( @@ -428,109 +537,34 @@ int cu_multMat( return err; } -extern "C" int cu_CmultMat( - cublasHandle_t hcublas, size_t m, size_t n, size_t k, +extern "C" int cu_Cmatmul( + void* hpblas, int32_t m, int32_t n, int32_t k, const cuComplex* alpha, const cuComplex* d_A, const cuComplex* d_B, const cuComplex* beta, cuComplex* d_C, int dagger ) { - return cu_multMat(hcublas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger); + return cu_matmul(hpblas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger); } -extern "C" int cu_ZmultMat( - cublasHandle_t hcublas, size_t m, size_t n, size_t k, +extern "C" int cu_Zmatmul( + void* hpblas, int32_t m, int32_t n, int32_t k, const cuDoubleComplex* alpha, const cuDoubleComplex* d_A, const cuDoubleComplex* d_B, const cuDoubleComplex* beta, cuDoubleComplex* d_C, int dagger ) { - return cu_multMat(hcublas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger); + return cu_matmul(hpblas, m, n, k, alpha, d_A, d_B, beta, d_C, dagger); } -template -int inverse( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, Number* d_A, - Number* d_Ainv, size_t n -) { - ENFORCE(hcusolver); - ENFORCE(d_A); - ENFORCE(d_Ainv); - - // compute buffer size and prep . memory - int lwork; - cusolverStatus_t cusolverStatus = - libnegf::cusolverDngetrf_bufferSize(hcusolver, n, n, d_A, n, &lwork); - ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); - - // prepare memory on the device - Number* d_LU; - cudaError_t cudaStatus = cudaMalloc((void**)&d_LU, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); - - int* d_pivot; - cudaStatus = cudaMalloc((void**)&d_pivot, n * sizeof(int)); - ENFORCE(cudaStatus == cudaSuccess); - int* d_info; - cudaStatus = cudaMalloc((void**)&d_info, sizeof(int)); - ENFORCE(cudaStatus == cudaSuccess); - // copy d_LU <- pdA - auto cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_A, 1, d_LU, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - - Number* d_work; - cudaStatus = cudaMalloc((void**)&d_work, lwork * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); - - // LU factorization of d_A , with partial pivoting and row - // interchanges ; row i is interchanged with row d_pivot ( i ); - cusolverStatus = libnegf::cusolverDngetrf( - hcusolver, n, n, d_LU, n, d_work, d_pivot, d_info - ); - ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); - LIBNEGF_CHECK_INFO(d_info); - - // use the LU factorization to solve the system d_LU * x = d_Ainv ; - // the solution overwrites d_Ainv - cusolverStatus = libnegf::cusolverDngetrs( - hcusolver, CUBLAS_OP_N, n, n, d_LU, n, d_pivot, d_Ainv, n, d_info - ); - ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); - LIBNEGF_CHECK_INFO(d_info); - - cudaStatus = cudaFree(d_pivot); - ENFORCE(cudaStatus == cudaSuccess); - cudaStatus = cudaFree(d_info); - ENFORCE(cudaStatus == cudaSuccess); - cudaStatus = cudaFree(d_work); - ENFORCE(cudaStatus == cudaSuccess); - cudaStatus = cudaFree(d_LU); - ENFORCE(cudaStatus == cudaSuccess); - - return cudaStatus; -} - -extern "C" int cu_Cinverse( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, cuComplex* d_A, - cuComplex* d_Ainv, size_t n -) { - return inverse(hcublas, hcusolver, d_A, d_Ainv, n); -} - -extern "C" int cu_Zinverse( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, cuDoubleComplex* d_A, - cuDoubleComplex* d_Ainv, size_t n -) { - return inverse(hcublas, hcusolver, d_A, d_Ainv, n); -} template int cu_kernelsum( Number* d_C, Number* alpha, Number* d_A, Number* beta, Number* d_B, - size_t size + size_t num_elements ) { - auto num_blocks = (size + BLOCK_SIZE - 1) / BLOCK_SIZE; + auto num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; - addKernel<<>>(d_C, *alpha, d_A, *beta, d_B, size); - ENFORCE(cudaGetLastError() == cudaSuccess); + addKernel<<>>(d_C, *alpha, d_A, *beta, d_B, num_elements); + CUDA_CHECK(cudaGetLastError()); return 0; } @@ -538,16 +572,16 @@ int cu_kernelsum( extern "C" int cu_Ckernelsum( cuComplex* d_C, cuComplex* alpha, cuComplex* d_A, cuComplex* beta, - cuComplex* d_B, size_t size + cuComplex* d_B, size_t num_elements ) { - return cu_kernelsum(d_C, alpha, d_A, beta, d_B, size); + return cu_kernelsum(d_C, alpha, d_A, beta, d_B, num_elements); } extern "C" int cu_Zkernelsum( cuDoubleComplex* d_C, cuDoubleComplex* alpha, cuDoubleComplex* d_A, - cuDoubleComplex* beta, cuDoubleComplex* d_B, size_t size + cuDoubleComplex* beta, cuDoubleComplex* d_B, size_t num_elements ) { - return cu_kernelsum(d_C, alpha, d_A, beta, d_B, size); + return cu_kernelsum(d_C, alpha, d_A, beta, d_B, num_elements); } @@ -556,9 +590,9 @@ extern "C" int cu_Zkernelsum( * or its complex-conjugate is used. * * The possible options are: - * * `dagger == 0`: compute C ≔ α A + β B - * * `dagger == 1`: compute C ≔ α A^* + β B - * * `dagger == 2`: compute C ≔ α A + B^* + * * `dagger == 0`: compute C = α A + β B + * * `dagger == 1`: compute C = α A^* + β B + * * `dagger == 2`: compute C = α A + B^* * * @param[in] m The number of rows of op(A) and op(B). * @param[in] n The number of columns of op(A) and op(B). @@ -566,7 +600,7 @@ extern "C" int cu_Zkernelsum( */ template int cu_matsum( - cublasHandle_t hcublas, size_t m, size_t n, const Number* alpha, + void *hpblas, int32_t m, int32_t n, const Number* alpha, const Number* d_A, const Number* beta, const Number* d_B, Number* d_C, int dagger ) { @@ -575,7 +609,7 @@ int cu_matsum( ENFORCE(d_C); ENFORCE(d_A != d_C || dagger == 0 || dagger == 2); ENFORCE(d_B != d_C || dagger == 1); - + cublasHandle_t hcublas = static_cast(hpblas); cublasStatus_t err; if(dagger == 0) { err = libnegf::cublasGeam( @@ -601,106 +635,120 @@ int cu_matsum( } extern "C" int cu_Cmatsum( - cublasHandle_t hcublas, size_t m, size_t n, const cuComplex* alpha, + void *hpblas, int32_t m, int32_t n, const cuComplex* alpha, const cuComplex* d_A, const cuComplex* beta, const cuComplex* d_B, cuComplex* d_C, int dagger ) { - return cu_matsum(hcublas, m, n, alpha, d_A, beta, d_B, d_C, dagger); + return cu_matsum(hpblas, m, n, alpha, d_A, beta, d_B, d_C, dagger); } extern "C" int cu_Zmatsum( - cublasHandle_t hcublas, size_t m, size_t n, const cuDoubleComplex* alpha, + void *hpblas, int32_t m, int32_t n, const cuDoubleComplex* alpha, const cuDoubleComplex* d_A, const cuDoubleComplex* beta, const cuDoubleComplex* d_B, cuDoubleComplex* d_C, int dagger ) { - return cu_matsum(hcublas, m, n, alpha, d_A, beta, d_B, d_C, dagger); + return cu_matsum(hpblas, m, n, alpha, d_A, beta, d_B, d_C, dagger); } -extern "C" int cu_Cinitmat(cuComplex* d_A, size_t nrow) { +extern "C" int cu_Cinitmat(cuComplex* d_A, int32_t nrow) { ENFORCE(d_A); - auto size = nrow * nrow; - auto num_blocks = (size + BLOCK_SIZE - 1) / BLOCK_SIZE; + size_t n64 = static_cast(nrow); + size_t num_elements = n64 * n64; + size_t num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; initKernel<<>>(d_A, nrow); - ENFORCE(cudaGetLastError() == cudaSuccess); + CUDA_CHECK(cudaGetLastError()); return 0; } -extern "C" int cu_Zinitmat(cuDoubleComplex* d_A, size_t nrow) { +extern "C" int cu_Zinitmat(cuDoubleComplex* d_A, int32_t nrow) { ENFORCE(d_A); - auto size = nrow * nrow; - auto num_blocks = (size + BLOCK_SIZE - 1) / BLOCK_SIZE; + size_t n64 = static_cast(nrow); + size_t num_elements = n64 * n64; + size_t num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; initKernel<<>>(d_A, nrow); - ENFORCE(cudaGetLastError() == cudaSuccess); + CUDA_CHECK(cudaGetLastError()); return 0; } + template::type> Real trace( - cublasHandle_t hcublas, Number* d_A, size_t nrow, void* h_mask, + Number* d_A, int32_t nrow, void* h_mask, int mask_present ) { - auto size = nrow * nrow; - auto num_blocks = (size + BLOCK_SIZE - 1) / BLOCK_SIZE; - Real* d_work; - cudaError_t cudaStatus = cudaMalloc((void**)&d_work, nrow * sizeof(Real)); - Real* d_iden; - cudaStatus = cudaMalloc((void**)&d_iden, nrow * sizeof(Real)); - ENFORCE(cudaStatus == cudaSuccess); + + size_t n64 = static_cast(nrow); + size_t num_elements = n64 * n64; + size_t num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; + + cudaError_t cudaStatus; + bool* d_mask; - cudaStatus = cudaMalloc((void**)&d_mask, nrow * sizeof(bool)); - ENFORCE(cudaStatus == cudaSuccess); - if(h_mask) { - cudaStatus = cudaMemcpy( - d_mask, h_mask, nrow * sizeof(bool), cudaMemcpyHostToDevice - ); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_mask, n64 * sizeof(bool)); + CUDA_CHECK(cudaStatus); + if (h_mask) { + cu_copyMatH2D(h_mask, d_mask, n64*sizeof(bool)); } - initArrayWithOnes<<>>(d_iden, nrow); - traceKernel<<>>( - d_A, nrow, d_work, d_mask, mask_present - ); + Real h_result; + Real *d_result; + Real *d_real_diag; + void *d_temp = nullptr; + size_t temp_bytes = 0; + cudaStatus = cudaMalloc(&d_result, sizeof(Real)); + cudaStatus = cudaMemset(static_cast(d_result), 0, sizeof(Real)); + CUDA_CHECK(cudaStatus); + cudaStatus = cudaMalloc(&d_real_diag, n64*sizeof(Real)); + CUDA_CHECK(cudaStatus); - Real result; - cublasStatus_t err = - libnegf::cublasDot(hcublas, nrow, d_iden, 1, d_work, 1, &result); - ENFORCE(err == CUBLAS_STATUS_SUCCESS); + getdiagKernel<<>>(d_A, nrow, d_real_diag, d_mask, mask_present); + + cub::DeviceReduce::Sum(nullptr, temp_bytes, d_real_diag, d_result, nrow); + + cudaStatus = cudaMalloc(&d_temp, temp_bytes); + CUDA_CHECK(cudaStatus); + + cub::DeviceReduce::Sum(d_temp, temp_bytes, d_real_diag, d_result, nrow); + + cudaStatus = cudaMemcpy(&h_result, d_result, sizeof(Real), cudaMemcpyDeviceToHost); - cudaStatus = cudaFree(d_work); - cudaStatus = cudaFree(d_iden); cudaStatus = cudaFree(d_mask); + cudaStatus = cudaFree(d_result); + cudaStatus = cudaFree(d_temp); + cudaStatus = cudaFree(d_real_diag); - return result; + return h_result; } extern "C" float cu_Ctrace( - cublasHandle_t hcublas, cuComplex* d_A, size_t nrow, void* h_mask, + void *hpblas, cuComplex* d_A, int32_t nrow, void* h_mask, int mask_present ) { - return trace(hcublas, d_A, nrow, h_mask, mask_present); + return trace(d_A, nrow, h_mask, mask_present); } extern "C" double cu_Ztrace( - cublasHandle_t hcublas, cuDoubleComplex* d_A, size_t nrow, void* h_mask, + void *hpblas, cuDoubleComplex* d_A, int32_t nrow, void* h_mask, int mask_present ) { - return trace(hcublas, d_A, nrow, h_mask, mask_present); + return trace(d_A, nrow, h_mask, mask_present); } extern "C" int cu_Cmatcopy( - cublasHandle_t hcublas, const cuComplex* d_A, cuComplex* d_B, + void *hpblas, const cuComplex* d_A, cuComplex* d_B, size_t num_elements ) { ENFORCE(d_A); ENFORCE(d_B); ENFORCE(num_elements <= INT_MAX); + cublasHandle_t hcublas = static_cast(hpblas); auto err = cublasCcopy(hcublas, num_elements, d_A, 1, d_B, 1); ENFORCE(err == CUBLAS_STATUS_SUCCESS); @@ -708,12 +756,13 @@ extern "C" int cu_Cmatcopy( } extern "C" int cu_Zmatcopy( - cublasHandle_t hcublas, cuDoubleComplex* d_A, cuDoubleComplex* d_B, + void *hpblas, cuDoubleComplex* d_A, cuDoubleComplex* d_B, size_t num_elements ) { ENFORCE(d_A); ENFORCE(d_B); ENFORCE(num_elements <= INT_MAX); + cublasHandle_t hcublas = static_cast(hpblas); auto err = cublasZcopy(hcublas, num_elements, d_A, 1, d_B, 1); ENFORCE(err == CUBLAS_STATUS_SUCCESS); @@ -721,123 +770,205 @@ extern "C" int cu_Zmatcopy( } extern "C" int -cu_Casum(cublasHandle_t hcublas, void* d_A, float* summ, size_t n) { +cu_Casum(void *hpblas, void* d_A, float* summ, int32_t n) { cuComplex* pdA = (cuComplex*)d_A; + cublasHandle_t hcublas = static_cast(hpblas); auto err = cublasScasum(hcublas, n, pdA, 1, summ); ENFORCE(err == CUBLAS_STATUS_SUCCESS); return err; } extern "C" int -cu_Zasum(cublasHandle_t hcublas, void* d_A, double* summ, size_t n) { +cu_Zasum(void *hpblas, void* d_A, double* summ, int32_t n) { cuDoubleComplex* pdA = (cuDoubleComplex*)d_A; + cublasHandle_t hcublas = static_cast(hpblas); auto err = cublasDzasum(hcublas, n, pdA, 1, summ); ENFORCE(err == CUBLAS_STATUS_SUCCESS); return err; } +/* ************************************************************************** + * Matrix inversion A^-1 + * ************************************************************************** + */ template +int inverse(void* hpsolver, Number* d_A, Number* d_Ainv, int32_t n +) { + ENFORCE(hpsolver); + ENFORCE(d_A); + ENFORCE(d_Ainv); + + size_t n64 = static_cast(n); + size_t num_elements = n64 * n64; + size_t bytecount = num_elements * sizeof(Number); + + cusolverDnHandle_t hcusolver = static_cast(hpsolver); + + cusolverStatus_t cusolverStatus; + cudaError_t cudaStatus; + + // compute buffer msize and prep . memory + int lwork; + cusolverStatus = + libnegf::cusolverDngetrf_bufferSize(hcusolver, n, n, d_A, n, &lwork); + ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); + + // prepare memory on the device + Number* d_LU; + cudaStatus = cudaMalloc((void**)&d_LU, bytecount); + CUDA_CHECK(cudaStatus); + + int* d_pivot; + cudaStatus = cudaMalloc((void**)&d_pivot, n64 * sizeof(int)); + CUDA_CHECK(cudaStatus); + int* d_info; + cudaStatus = cudaMalloc((void**)&d_info, sizeof(int)); + CUDA_CHECK(cudaStatus); + // copy d_LU <- pdA + cu_copyMatD2D(d_LU, d_A, bytecount); + + Number* d_work; + size_t lwork_size = static_cast(lwork) * sizeof(Number); + cudaStatus = cudaMalloc((void**)&d_work, lwork_size); + CUDA_CHECK(cudaStatus); + + // LU factorization of d_A , with partial pivoting and row + // interchanges ; row i is interchanged with row d_pivot ( i ); + cusolverStatus = libnegf::cusolverDngetrf( + hcusolver, n, n, d_LU, n, d_work, d_pivot, d_info + ); + ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); + LIBNEGF_CHECK_INFO(d_info); + + // use the LU factorization to solve the system d_LU * x = d_Ainv ; + // the solution overwrites d_Ainv + cusolverStatus = libnegf::cusolverDngetrs( + hcusolver, CUBLAS_OP_N, n, n, d_LU, n, d_pivot, d_Ainv, n, d_info + ); + ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); + LIBNEGF_CHECK_INFO(d_info); + + cudaStatus = cudaFree(d_pivot); + CUDA_CHECK(cudaStatus); + cudaStatus = cudaFree(d_info); + CUDA_CHECK(cudaStatus); + cudaStatus = cudaFree(d_work); + CUDA_CHECK(cudaStatus); + cudaStatus = cudaFree(d_LU); + CUDA_CHECK(cudaStatus); + + return cudaStatus; +} + +extern "C" int cu_Cinverse( + void* hpsolver, cuComplex* d_A, + cuComplex* d_Ainv, int32_t n +) { + return inverse(hpsolver, d_A, d_Ainv, n); +} + +extern "C" int cu_Zinverse( + void * hpsolver, cuDoubleComplex* d_A, + cuDoubleComplex* d_Ainv, int32_t n +) { + return inverse(hpsolver, d_A, d_Ainv, n); +} + +/* *************************************************************** + * Sanchez-Sanchez algorithm + * **************************************************************** + */ +template::type> int decimation( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, Number* h_Go_out, - Number* h_Ao_in, Number* h_Bo_in, Number* h_Co_in, size_t n, int tf32, - int* ncyc, float SGFACC + void* hpblas, void* hpsolver, Number* h_Go_out, + Number* h_Ao_in, Number* h_Bo_in, Number* h_Co_in, int32_t n, + int* ncyc, Real accuracy, bool log ) { ENFORCE(h_Go_out); ENFORCE(h_Ao_in); ENFORCE(h_Bo_in); ENFORCE(h_Co_in); - ENFORCE(tf32 == 0 || tf32 == 1); ENFORCE(ncyc); - ENFORCE(SGFACC > 0.0); + ENFORCE(accuracy > 0.0); - auto num_elements = n * n; - auto num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; + cublasHandle_t hcublas = static_cast(hpblas); + cusolverDnHandle_t hcusolver = static_cast(hpsolver); + + cudaError_t cudaStatus; + cublasStatus_t err; + + size_t n64 = static_cast(n); + size_t num_elements = n64 * n64; + size_t bytecount = num_elements * sizeof(Number); + size_t num_blocks = (num_elements + BLOCK_SIZE - 1) / BLOCK_SIZE; - using Real = typename get_real::type; auto one = Number{Real{1}}; auto mone = Number{-Real{1}}; auto zero = Number{Real{0}}; Number* d_Ao; - cudaError_t cudaStatus = - cudaMalloc((void**)&d_Ao, num_elements * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Ao, bytecount); + CUDA_CHECK(cudaStatus); Number* d_Bo; - cudaStatus = cudaMalloc((void**)&d_Bo, num_elements * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Bo, bytecount); + CUDA_CHECK(cudaStatus); Number* d_Co; - cudaStatus = cudaMalloc((void**)&d_Co, num_elements * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); - - cudaStatus = cudaMemcpy( - d_Ao, h_Ao_in, n * n * sizeof(Number), cudaMemcpyHostToDevice - ); - ENFORCE(cudaStatus == cudaSuccess); - cudaStatus = cudaMemcpy( - d_Bo, h_Bo_in, n * n * sizeof(Number), cudaMemcpyHostToDevice - ); - ENFORCE(cudaStatus == cudaSuccess); - cudaStatus = cudaMemcpy( - d_Co, h_Co_in, n * n * sizeof(Number), cudaMemcpyHostToDevice - ); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Co, bytecount); + CUDA_CHECK(cudaStatus); - cublasStatus_t cublasStatus = - cublasSetPointerMode(hcublas, CUBLAS_POINTER_MODE_HOST); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - - if(tf32 == 1) { - cublasStatus = cublasSetMathMode(hcublas, CUBLAS_TENSOR_OP_MATH); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - } + cu_copyMatH2D(h_Ao_in, d_Ao, bytecount); + cu_copyMatH2D(h_Bo_in, d_Bo, bytecount); + cu_copyMatH2D(h_Co_in, d_Co, bytecount); Number* d_Ao_s; - cudaStatus = cudaMalloc((void**)&d_Ao_s, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Ao_s, bytecount); + CUDA_CHECK(cudaStatus); Number* d_C1; - cudaStatus = cudaMalloc((void**)&d_C1, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_C1, bytecount); + CUDA_CHECK(cudaStatus); Number* d_Go; - cudaStatus = cudaMalloc((void**)&d_Go, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Go, bytecount); + CUDA_CHECK(cudaStatus); int* d_pivot; - cudaStatus = cudaMalloc((void**)&d_pivot, n * sizeof(int)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_pivot, n64 * sizeof(int)); + CUDA_CHECK(cudaStatus); Number* d_T; - cudaStatus = cudaMalloc((void**)&d_T, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_T, bytecount); + CUDA_CHECK(cudaStatus); Number* d_Self; - cudaStatus = cudaMalloc((void**)&d_Self, n * n * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + cudaStatus = cudaMalloc((void**)&d_Self, bytecount); + CUDA_CHECK(cudaStatus); int* d_info; cudaStatus = cudaMalloc((void**)&d_info, sizeof(int)); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); + Real* d_norm; + cudaStatus = cudaMalloc((void**)&d_norm, sizeof(Real)); + CUDA_CHECK(cudaStatus); int lwork; cusolverStatus_t cusolverStatus = - libnegf::cusolverDngetrf_bufferSize(hcusolver, n, n, d_Self, n, &lwork); + libnegf::cusolverDngetrf_bufferSize(hcusolver, n, n, d_Self, n, &lwork); ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); Number* d_work; cudaStatus = cudaMalloc((void**)&d_work, lwork * sizeof(Number)); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); - cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_Ao, 1, d_Ao_s, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + cu_copyMatD2D(d_Ao_s, d_Ao, bytecount); constexpr auto MAX_NUM_ITERATIONS = 300; bool okCo = false; for(int i1 = 1; i1 <= MAX_NUM_ITERATIONS; i1++) { *ncyc = i1; + // Compute Go = Ao^-1 initKernel<<>>(d_Go, n); - cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_Ao, 1, d_Self, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + cu_copyMatD2D(d_Self, d_Ao, bytecount); cusolverStatus = libnegf::cusolverDngetrf( hcusolver, n, n, d_Self, n, d_work, d_pivot, d_info @@ -850,24 +981,31 @@ int decimation( ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); LIBNEGF_CHECK_INFO(d_info); - cublasStatus = libnegf::cublasGemm( + // compute C1 = Co Go Co + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &one, d_Go, n, d_Co, n, &zero, d_T, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); - cublasStatus = libnegf::cublasGemm( + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &one, d_Co, n, d_T, n, &zero, d_C1, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); + + // convergence check on |C1| < tol + cudaStatus = cudaMemset(d_norm, 0, sizeof(Real)); + CUDA_CHECK(cudaStatus); - Real summ; - cublasStatus = libnegf::cublasAsum(hcublas, n * n, d_C1, 1, &summ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - // printf("loop it= %d , summ= %f \n ", i1, summ); + maxabs_kernel<<>> + (d_C1, n * n, d_norm); + CUDA_CHECK(cudaGetLastError()); + Real h_norm; + cudaStatus = cudaMemcpy(&h_norm, d_norm, sizeof(Real), cudaMemcpyDeviceToHost); + if (log){printf("loop it= %d , maxval|C1|= %g \n ", i1, h_norm);} - if(summ <= SGFACC) { + if(h_norm <= accuracy) { if(okCo) { break; } else { @@ -877,48 +1015,53 @@ int decimation( okCo = false; } - cublasStatus = libnegf::cublasGemm( + // compute Self = Bo Go Co + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &one, d_Bo, n, d_T, n, &zero, d_Self, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); - cublasStatus = - libnegf::cublasAxpy(hcublas, n * n, &mone, d_Self, 1, d_Ao_s, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - cublasStatus = - libnegf::cublasAxpy(hcublas, n * n, &mone, d_Self, 1, d_Ao, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + // compute Ao_s = Ao_s - Self + err = libnegf::cublasAxpy(hcublas, n * n, &mone, d_Self, 1, d_Ao_s, 1); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); - cublasStatus = libnegf::cublasGemm( + // compute Ao = Ao - Self + err = libnegf::cublasAxpy(hcublas, n * n, &mone, d_Self, 1, d_Ao, 1); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); + + // compute Ao = Ao - Co Go Bo + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &one, d_Go, n, d_Bo, n, &zero, d_T, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); - cublasStatus = libnegf::cublasGemm( + ENFORCE(err == CUBLAS_STATUS_SUCCESS); + + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &mone, d_Co, n, d_T, n, &one, d_Ao, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); - cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_C1, 1, d_Co, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + // Co = C1 + cu_copyMatD2D(d_Co, d_C1, bytecount); - cublasStatus = libnegf::cublasGemm( + // compute C1 = Bo Go Bo + err = libnegf::cublasGemm( hcublas, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n, &one, d_Bo, n, d_T, n, &zero, d_C1, n ); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + ENFORCE(err == CUBLAS_STATUS_SUCCESS); - cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_C1, 1, d_Bo, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + // compute Bo = C1 + cu_copyMatD2D(d_Bo, d_C1, bytecount); } ENFORCE(*ncyc < MAX_NUM_ITERATIONS); + // compute Go = Ao_s^-1 initKernel<<>>(d_Go, n); - cublasStatus = libnegf::cublasCopy(hcublas, n * n, d_Ao_s, 1, d_Self, 1); - ENFORCE(cublasStatus == CUBLAS_STATUS_SUCCESS); + cu_copyMatD2D(d_Self, d_Ao_s, bytecount); cusolverStatus = libnegf::cusolverDngetrf( hcusolver, n, n, d_Self, n, d_work, d_pivot, d_info ); @@ -930,67 +1073,67 @@ int decimation( ENFORCE(cusolverStatus == CUSOLVER_STATUS_SUCCESS); LIBNEGF_CHECK_INFO(d_info); - cudaStatus = cudaMemcpy( - h_Go_out, d_Go, n * n * sizeof(Number), cudaMemcpyDeviceToHost - ); - ENFORCE(cudaStatus == cudaSuccess); + // copy to output to host + cu_copyMatD2H(h_Go_out, d_Go, bytecount); cudaStatus = cudaFree(d_pivot); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_info); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Ao); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Bo); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Co); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Go); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Ao_s); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_C1); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_T); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_Self); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); cudaStatus = cudaFree(d_work); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); + cudaStatus = cudaFree(d_norm); + CUDA_CHECK(cudaStatus); return cudaStatus; } extern "C" int cu_Cdecimation( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, cuComplex* h_Go_out, - cuComplex* h_Ao_in, cuComplex* h_Bo_in, cuComplex* h_Co_in, size_t n, - int tf32, int* ncyc, float SGFACC + void* hcublas, void* hcusolver, cuComplex* h_Go_out, + cuComplex* h_Ao_in, cuComplex* h_Bo_in, cuComplex* h_Co_in, int32_t n, + int* ncyc, float accuracy, bool log = false ) { return decimation( - hcublas, hcusolver, h_Go_out, h_Ao_in, h_Bo_in, h_Co_in, n, tf32, ncyc, - SGFACC + hcublas, hcusolver, h_Go_out, h_Ao_in, h_Bo_in, h_Co_in, n, ncyc, + accuracy, log ); } extern "C" int cu_Zdecimation( - cublasHandle_t hcublas, cusolverDnHandle_t hcusolver, + void* hcublas, void* hcusolver, cuDoubleComplex* h_Go_out, cuDoubleComplex* h_Ao_in, - cuDoubleComplex* h_Bo_in, cuDoubleComplex* h_Co_in, size_t n, int tf32, - int* ncyc, double SGFACC + cuDoubleComplex* h_Bo_in, cuDoubleComplex* h_Co_in, int32_t n, + int* ncyc, double accuracy, bool log = false ) { return decimation( - hcublas, hcusolver, h_Go_out, h_Ao_in, h_Bo_in, h_Co_in, n, tf32, ncyc, - SGFACC + hcublas, hcusolver, h_Go_out, h_Ao_in, h_Bo_in, h_Co_in, n, ncyc, + accuracy, log ); } extern "C" int cu_meminfo(size_t* freemem, size_t* totalmem) { //auto cudaStatus = cudaDeviceSynchronize(); - //printf("STATUS: %d\n",cudaStatus); - //ENFORCE(cudaStatus == cudaSuccess); + //printf("STATUS: %d\n",cudaStatus); + //CUDA_CHECK(cudaStatus); auto cudaStatus = cudaMemGetInfo(freemem, totalmem); - ENFORCE(cudaStatus == cudaSuccess); + CUDA_CHECK(cudaStatus); return cudaStatus; } diff --git a/src/iterative.F90 b/src/iterative.F90 index 0eeb5997..506af787 100644 --- a/src/iterative.F90 +++ b/src/iterative.F90 @@ -1741,6 +1741,9 @@ subroutine calculate_transmissions(negf, Ec, SelfEneR, tun_proj, tun_mat) if (nt.gt.1) then call calculate_Gr_tridiag_blocks(negf,ESH,gsmr,Gr,2,nt) end if +#:if defined("GPU") + call waitForGPU() +#:endif else ! When more contacts are present sometimes we can re-use previous GF ! if nt1 > nt extend the Gr calculation @@ -1749,6 +1752,9 @@ subroutine calculate_transmissions(negf, Ec, SelfEneR, tun_proj, tun_mat) nt = nt1 endif end if +#:if defined("GPU") + call waitForGPU() +#:endif call calculate_single_transmission_N_contacts(negf,nit,nft,ESH,SelfEneR,cblk,negf%tun_proj,gsmr,Gr,tun) @@ -1846,9 +1852,6 @@ subroutine calculate_transmissions_and_dos(negf, Ec, SelfEneR, GS, tun_proj, tun call calculate_Gr_tridiag_blocks(negf,ESH,gsmr,Gr,1) call calculate_Gr_tridiag_blocks(negf,ESH,gsmr,Gr,2,nbl) -#:if defined("GPU") - call waitForGPU() -#:endif !Computation of transmission(s) between contacts ni(:) -> nf(:) #:if defined("GPU") call waitForGPU() diff --git a/src/iterative_cpu.F90 b/src/iterative_cpu.F90 index 8b2a355b..662cda9e 100644 --- a/src/iterative_cpu.F90 +++ b/src/iterative_cpu.F90 @@ -255,7 +255,7 @@ subroutine calculate_Gn_tridiag_blocks(negf,ESH,SelfEneR,frm,ref,struct,gsmr,Gr, ! Add contact self-energies do j=1,ncont frmdiff = frm(j) - frm(ref) - if (j.NE.ref .AND. ABS(frmdiff).GT.EPS) THEN + if (j.NE.ref .AND. ABS(frmdiff).GT.EPS10) THEN cb=struct%cblk(j) ! block corresponding to contact j call zspectral(SelfEneR(j),SelfEneR(j),0,Gam) Sigma_n(cb,cb)%val = Sigma_n(cb,cb)%val + frmdiff*Gam%val diff --git a/src/iterative_gpu.F90 b/src/iterative_gpu.F90 index da075de1..f0c0ad57 100644 --- a/src/iterative_gpu.F90 +++ b/src/iterative_gpu.F90 @@ -51,98 +51,98 @@ module iterative_gpu interface calculate_gsmr_blocks #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure calculate_gsmr_blocks_${KIND}$ #:endfor end interface calculate_gsmr_blocks interface calculate_Gr_tridiag_blocks #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure calculate_Gr_tridiag_blocks_${KIND}$ #:endfor end interface calculate_Gr_tridiag_blocks interface calculate_Gn_tridiag_blocks #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure calculate_Gn_tridiag_blocks_${KIND}$ #:endfor end interface calculate_Gn_tridiag_blocks interface calculate_single_transmission_2_contacts #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure calculate_single_transmission_2_contacts_${KIND}$ #:endfor end interface calculate_single_transmission_2_contacts interface calculate_single_transmission_N_contacts #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure calculate_single_transmission_N_contacts_${KIND}$ #:endfor end interface calculate_single_transmission_N_contacts interface get_tun_mask #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure get_tun_mask_${KIND}$ #:endfor end interface get_tun_mask interface check_convergence_trid #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure check_convergence_trid_${KIND}$ #:endfor end interface check_convergence_trid interface check_convergence_vec #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure check_convergence_vec_${KIND}$ #:endfor end interface check_convergence_vec - interface add_sigma_n + interface add_sigma_n #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure add_sigma_n_${KIND}$ #:endfor end interface add_sigma_n - interface init_tridiag_blk + interface init_tridiag_blk #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure init_tridiag_blk_${KIND}$ #:endfor end interface init_tridiag_blk - interface destroy_tridiag_blk + interface destroy_tridiag_blk #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure destroy_tridiag_blk_${KIND}$ #:endfor end interface destroy_tridiag_blk interface allocate_blk_dns #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure allocate_blk_dns_${KIND}$ #:endfor end interface allocate_blk_dns interface deallocate_blk_dns #:for PREC in PRECISIONS - #:set KIND = PREC_ABBREVS[PREC] + #:set KIND = PREC_ABBREVS[PREC] module procedure deallocate_blk_dns_${KIND}$ #:endfor end interface deallocate_blk_dns contains - + #:def calculate_gsmr_blocks_template(KIND, CTYPE, MTYPE, CUDATYPE) subroutine calculate_gsmr_blocks_${KIND}$(negf,ESH,sbl,ebl,gsmr,keep_gsmr) @@ -158,9 +158,9 @@ subroutine calculate_gsmr_blocks_${KIND}$(negf,ESH,sbl,ebl,gsmr,keep_gsmr) !Work type(CublasHandle) :: hh type(CusolverDnHandle) :: hhsol - complex(${KIND}$), parameter :: one = cmplx(1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: mone = cmplx(-1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: zero = cmplx(0.0, 0.0, ${KIND}$) + complex(${KIND}$) :: one = cmplx(1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: mone = cmplx(-1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: zero = cmplx(0.0, 0.0, ${KIND}$) type(${MTYPE}$) :: work1, work2 integer :: nrow integer :: i, nbl, istat @@ -239,9 +239,9 @@ subroutine calculate_Gr_tridiag_blocks_${KIND}$(negf,ESH,gsmr,Gr,sbl,ebl) type(CusolverDnHandle) ::hhsol integer :: i,nbl type(${MTYPE}$), target :: work1, work2, work3 - complex(${KIND}$), parameter :: one = cmplx(1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: mone = cmplx(-1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: zero = cmplx(0.0, 0.0, ${KIND}$) + complex(${KIND}$) :: one = cmplx(1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: mone = cmplx(-1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: zero = cmplx(0.0, 0.0, ${KIND}$) integer :: istat nbl = size(ESH,1) @@ -433,9 +433,9 @@ subroutine calculate_Gn_tridiag_blocks_${KIND}$(negf,ESH,SelfEneR,frm,ref,struct !Work type(CublasHandle) :: hh - complex(${KIND}$), parameter :: one = cmplx(1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: minusone = cmplx(-1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: zero = cmplx(0.0, 0.0, ${KIND}$) + complex(${KIND}$) :: one = cmplx(1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: minusone = cmplx(-1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: zero = cmplx(0.0, 0.0, ${KIND}$) type(${MTYPE}$), dimension(:,:), allocatable :: Sigma_n type(${MTYPE}$) :: work1, Gam complex(${KIND}$) :: frmdiff @@ -445,7 +445,6 @@ subroutine calculate_Gn_tridiag_blocks_${KIND}$(negf,ESH,SelfEneR,frm,ref,struct ncont = struct%num_conts nbl = struct%num_PLs hh = negf%hcublas - !build Sigma_n from SelfEneR call allocate_blk_dns(Sigma_n, nbl) call init_tridiag_blk(Sigma_n, ESH) @@ -453,22 +452,21 @@ subroutine calculate_Gn_tridiag_blocks_${KIND}$(negf,ESH,SelfEneR,frm,ref,struct ! Add interaction self-energies call add_sigma_n(negf, Sigma_n) - call copy_trid_toGPU(Sigma_n) ! Add contact self-energies do j=1,ncont frmdiff = cmplx(frm(j) - frm(ref), 0.0, ${KIND}$) - if (j.NE.ref .AND. ABS(frmdiff).GT.EPS) THEN + if (j.ne.ref .and. abs(frmdiff).gt.EPS10) then cb=struct%cblk(j) ! block corresponding to contact j call zspectral(SelfEneR(j),SelfEneR(j),0,Gam) - call createGPU(Gam) - call copyToGPU(Gam) - call matsum_gpu(hh, one, Sigma_n(cb,cb), frmdiff, Gam, Sigma_n(cb,cb)) - call destroyAll(Gam) + Sigma_n(cb,cb)%val = Sigma_n(cb,cb)%val + frmdiff * Gam%val + call destroy(Gam) endif end do - call calculate_sigma_n() + call copy_trid_toGPU(Sigma_n) + call calculate_sigma_n() + ! Gr * Sigma_n * Gr^+ call createAll(work1,Sigma_n(1,1)%nrow, Gr(1,1)%nrow) call matmul_gpu(hh, one, Sigma_n(1,1), Gr(1,1), zero, work1, 'dag_2nd') call matmul_gpu(hh, one, Gr(1,1), work1, zero, Gn(1,1)) @@ -607,7 +605,7 @@ subroutine add_sigma_n_${KIND}$(negf, sigma_n) type(TInteractionNode), pointer :: it it => negf%interactList%first - + do while (associated(it)) call it%inter%add_sigma_n(sigma_n, negf%iEloc, negf%iKloc, negf%spin) it => it%next @@ -632,11 +630,11 @@ subroutine calculate_single_transmission_2_contacts_${KIND}$(negf,ni,nf,ESH,Self type(CublasHandle) :: hh Integer :: ct1, bl1 Type(${MTYPE}$) :: work1, work2, GAM1_dns, TRS, AA - complex(${KIND}$), parameter :: j = cmplx(0.0, 1.0, ${KIND}$) ! CMPX unity - complex(${KIND}$), parameter :: mj = cmplx(0.0, -1.0, ${KIND}$) - complex(${KIND}$), parameter :: one = cmplx(1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: mone = cmplx(-1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: zero = cmplx(0.0, 0.0, ${KIND}$) + complex(${KIND}$) :: j = cmplx(0.0, 1.0, ${KIND}$) ! CMPX unity + complex(${KIND}$) :: mj = cmplx(0.0, -1.0, ${KIND}$) + complex(${KIND}$) :: one = cmplx(1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: mone = cmplx(-1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: zero = cmplx(0.0, 0.0, ${KIND}$) if (size(cblk).gt.2) then write(*,*) "ERROR: calculate_single_transmission_2_contacts is valid only for 2 contacts" @@ -714,9 +712,9 @@ subroutine calculate_single_transmission_N_contacts_${KIND}$(negf,ni,nf,ESH,Self Integer :: ct1, ct2, bl1, bl2, i, nbl Type(${MTYPE}$) :: work1, work2, GAM1_dns, GAM2_dns, TRS Real(${KIND}$) :: max - complex(${KIND}$), parameter :: one = cmplx(1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: mone = cmplx(-1.0, 0.0, ${KIND}$) - complex(${KIND}$), parameter :: zero = cmplx(0.0, 0.0, ${KIND}$) + complex(${KIND}$) :: one = cmplx(1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: mone = cmplx(-1.0, 0.0, ${KIND}$) + complex(${KIND}$) :: zero = cmplx(0.0, 0.0, ${KIND}$) !Arrange contacts in way that order between first and second is always the !same (always ct1 < ct2) @@ -891,7 +889,7 @@ subroutine destroy_tridiag_blk_${KIND}$(M, str) nbl=size(M,1) do i=1,nbl - if (allocated(M(i,i)%val)) then + if (allocated(M(i,i)%val)) then call destroyAll(M(i,i)) end if end do diff --git a/src/lib_param.F90 b/src/lib_param.F90 index a01cec95..6a448ccb 100644 --- a/src/lib_param.F90 +++ b/src/lib_param.F90 @@ -22,7 +22,7 @@ module lib_param use iso_c_binding - use ln_precision, only : dp + use ln_precision, only : dp, EPS15, set_drop use globals use mat_def use ln_structure, only : TStruct_info, TBasisCenters, TNeighbourMap @@ -88,11 +88,11 @@ module lib_param #:if defined("GPU") !Type definition of cublas and cusolver handles (no module defs?) type, bind(C) :: cublasHandle - type(c_ptr) :: handle + type(c_ptr) :: ptr end type cublasHandle type, bind(C) :: cusolverDnHandle - type(c_ptr) :: handle + type(c_ptr) :: ptr end type cusolverDnHandle public :: cusolverDnHandle @@ -666,6 +666,8 @@ subroutine set_defaults(negf) call init_cache_space(negf, 'G_r') call init_cache_space(negf, 'G_n') + call set_drop(EPS15) + end subroutine set_defaults diff --git a/src/sparskit/CMakeLists.txt b/src/sparskit/CMakeLists.txt index e69de29b..8b137891 100644 --- a/src/sparskit/CMakeLists.txt +++ b/src/sparskit/CMakeLists.txt @@ -0,0 +1 @@ + diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f7da815c..489a3d45 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,127 +1,170 @@ +# On some HPC systems (e.g., LUMI), front-end nodes do not feature GPUs or the front-end nodes may +# have an instruction set architecture different from compute nodes (e.g., on Fugaku with its Intel +# front-end and Arm64FX compute nodes). On these systems, GPU tests or all tests, respectively, have +# to be submitted with mpiexec/mpirun. +option(LIBNEGF_TEST_WITH_MPIEXEC "Always launch tests with mpiexec/mpirun" OFF) + # Test input files may be large and are therefore managed by git-lfs. find_program(GIT_LFS git-lfs) if(NOT GIT_LFS) - message(FATAL_ERROR - "Tests cannot be run because git-lfs is missing. " - "Disable testing or install git-lfs." + message(FATAL_ERROR "Tests cannot be run because git-lfs is missing. " + "Disable testing or install git-lfs." ) endif() - # List of active tets. Each test indicates a subdirectory. set(test-directories - f90init - c_mpi_init - cpp_mpi_init - c_int - c_int_elph_deph - c_int_file - f90int - f90int_file - f90_transmission - f90elph_deph - f90read_hs - f90Si2x2 - f90Si_nin - f90Si_nin_40pl - ) + f90init + c_mpi_init + cpp_mpi_init + c_int + c_int_elph_deph + c_int_file + f90int + f90int_file + f90_transmission + f90elph_deph + f90read_hs + f90Si2x2 + f90Si_nin + f90Si_nin_40pl +) if(WITH_TRANSPORT_GPU) list(APPEND test-directories testCUDA) list(APPEND test-directories testCUDA_decimation) endif() -# Define a function wich copies over the content of the tests. -# This is used to support out-of-source build, as the test -# directory contain input files which we don't want to specify. +# Define a function wich copies over the content of the tests. This is used to support out-of-source +# build, as the test directory contain input files which we don't want to specify. function(transfer_test_data testname) - if("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}") - add_custom_target(${testname}-data-transfer) - else() - add_custom_target( - ${testname}-data-transfer - COMMENT "Copying ${CMAKE_CURRENT_SOURCE_DIR} for out-of-source build." - COMMAND cp -f ${CMAKE_CURRENT_SOURCE_DIR}/* ${CMAKE_CURRENT_BINARY_DIR}/) - endif() - add_dependencies(${testname} ${testname}-data-transfer) + if("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}") + add_custom_target(${testname}-data-transfer) + else() + add_custom_target( + ${testname}-data-transfer + COMMENT "Copying ${CMAKE_CURRENT_SOURCE_DIR} for out-of-source build." + COMMAND cp -f ${CMAKE_CURRENT_SOURCE_DIR}/* ${CMAKE_CURRENT_BINARY_DIR}/ + ) + endif() + add_dependencies(${testname} ${testname}-data-transfer) endfunction() function(unzip_test_data testname) - find_program(XZ NAMES xz) - if(NOT XZ) - message(FATAL_ERROR "cannot run test '${testname}' because xz executable is missing") - endif() - - if("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}") - add_custom_target(${testname}-data-unzip) - else() - add_custom_target( - ${testname}-data-unzip - COMMENT "unzip ${CMAKE_CURRENT_SOURCE_DIR} for out-of-source build." - COMMAND tar -xJf ${CMAKE_CURRENT_SOURCE_DIR}/*.tar.xz -C ${CMAKE_CURRENT_BINARY_DIR}) - endif() - add_dependencies(${testname}-data-unzip ${testname}-data-transfer) - add_dependencies(${testname} ${testname}-data-unzip) -endfunction() + find_program(XZ NAMES xz) + if(NOT XZ) + message(FATAL_ERROR "cannot run test '${testname}' because xz executable is missing") + endif() -# A function to setup C++ tests. -function(setup_cpp_test testname label source) - message(STATUS "add test: " ${testname} " " ${source} " " ${label}) - add_executable(${testname} ${source}) - transfer_test_data(${testname}) - target_link_libraries(${testname} MPI::MPI_CXX) - target_link_libraries(${testname} negf) - target_link_libraries(${testname} LAPACK::LAPACK) - target_include_directories(${testname} PRIVATE ${CMAKE_BINARY_DIR}/src/api) - add_test( - NAME ${testname} - COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ./${testname}) - set_property(TEST ${testname} PROPERTY LABELS ${label}) + if("${CMAKE_SOURCE_DIR}" STREQUAL "${CMAKE_BINARY_DIR}") + add_custom_target(${testname}-data-unzip) + else() + add_custom_target( + ${testname}-data-unzip + COMMENT "unzip ${CMAKE_CURRENT_SOURCE_DIR} for out-of-source build." + COMMAND tar -xJf ${CMAKE_CURRENT_SOURCE_DIR}/*.tar.xz -C ${CMAKE_CURRENT_BINARY_DIR} + ) + endif() + add_dependencies(${testname}-data-unzip ${testname}-data-transfer) + add_dependencies(${testname} ${testname}-data-unzip) endfunction() -# A function to setup C tests. -function(setup_c_test testname label source) - message(STATUS "add test: " ${testname} " " ${source} " " ${label}) - add_executable(${testname} ${source} ${ARGN}) - transfer_test_data(${testname}) +function(setup_test testname source_main) + set(options GPUONLY MPI) + set(multi_value_args SOURCES LABELS) + cmake_parse_arguments(PARSE_ARGV 2 LIBNEGF_TEST "${options}" "" "${multi_value_args}") + + if(LIBNEGF_TEST_UNPARSED_ARGUMENTS) + message(WARNING "could not parse the following arguments: ${LIBNEGF_TEST_UNPARSED_ARGUMENTS}") + endif() + if(LIBNEGF_TEST_KEYWORDS_MISSING_VALUES) + message( + WARNING + "the following keyword arguments are missing values: ${LIBNEGF_TEST_KEYWORDS_MISSING_VALUES}" + ) + endif() + + list(PREPEND LIBNEGF_TEST_SOURCES ${source_main}) + if(LIBNEGF_TEST_GPUONLY) + list(APPEND LIBNEGF_TEST_LABELS gpu-only) + endif() + if(LIBNEGF_TEST_MPI) + list(APPEND LIBNEGF_TEST_LABELS mpi) + endif() + + # detect implementation language + if(source_main MATCHES "[.]c$") + set(LIBNEGF_TEST_LANG C) + elseif(source_main MATCHES "[.]cpp$") + set(LIBNEGF_TEST_LANG CXX) + elseif(source_main MATCHES "[.](f|F)90$") + set(LIBNEGF_TEST_LANG Fortran) + else() + message(SEND_ERROR "unknown file extension of test source ${source_main}") + endif() + + message(DEBUG "add test: ${testname} ${LIBNEGF_TEST_SOURCES} ${LIBNEGF_TEST_LABELS}") + + # set up test compilation + add_executable(${testname} ${LIBNEGF_TEST_SOURCES}) + + if(LIBNEGF_TEST_LANG STREQUAL C) target_link_libraries(${testname} MPI::MPI_C) - target_link_libraries(${testname} negf) - target_link_libraries(${testname} LAPACK::LAPACK) target_include_directories(${testname} PRIVATE ${CMAKE_BINARY_DIR}/src/api) - add_test(${testname} ${testname}) - set_property(TEST ${testname} PROPERTY LABELS ${label}) -endfunction() - -# A function to setup F90 tests. -function(setup_f90_test testname label source) - message(STATUS "add test: " ${testname} " " ${source} " " ${label}) - add_executable(${testname} ${source} ${ARGN}) - transfer_test_data(${testname}) + elseif(LIBNEGF_TEST_LANG STREQUAL CXX) + target_link_libraries(${testname} MPI::MPI_CXX) + target_include_directories(${testname} PRIVATE ${CMAKE_BINARY_DIR}/src/api) + elseif(LIBNEGF_TEST_LANG STREQUAL Fortran) target_link_libraries(${testname} MPI::MPI_Fortran) - target_link_libraries(${testname} negf) - target_link_libraries(${testname} LAPACK::LAPACK) target_include_directories(${testname} PRIVATE ${BUILD_MOD_DIR}) - add_test(${testname} ${testname}) - set_property(TEST ${testname} PROPERTY LABELS ${label}) -endfunction() + endif() + target_link_libraries(${testname} negf) + target_link_libraries(${testname} LAPACK::LAPACK) -# A function to setup F90 tests with MPI. -function(setup_f90_mpi_test testname label source) - message(STATUS "add test: " ${testname} " " ${source} " " ${label}) - add_executable(${testname} ${source} ${ARGN}) - transfer_test_data(${testname}) - target_link_libraries(${testname} MPI::MPI_Fortran) - target_link_libraries(${testname} negf) - target_link_libraries(${testname} LAPACK::LAPACK) - target_include_directories(${testname} PRIVATE ${BUILD_MOD_DIR}) - add_test( - NAME ${testname} - COMMAND env OMP_NUM_THREADS=1 ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ./${testname}) - set_property(TEST ${testname} PROPERTY LABELS ${label}) -endfunction() + # determine how to execute the test + if(${LIBNEGF_TEST_MPI}) + set(num_procs ${TEST_NUM_PROCS}) + set(cmd + # When running on a user's PC, it is possible to thrash the machine by having one MPI task + # per core and one OpenMP thread per task and per core (the default). The restriction below + # aims to avoid this scenario. + # + # `env` is wrapping the call to mpiexec/mpirun to avoid the job name `env` when + # ${MPIEXEC_EXECUTABLE} is a batch scheduler. This works with Slurm but some batch + # schedulers (e.g., OAR) filter the set of environment variables (the OAR documentation [1] + # even mentions `OMP_NUM_THREADS` explicitly). + # + # [1] https://oar.imag.fr/wiki:passing_environment_variables_to_openmpi_nodes + env + OMP_NUM_THREADS=${TEST_NUM_THREADS} + ${MPIEXEC_EXECUTABLE} + ${MPIEXEC_NUMPROC_FLAG} + ${num_procs} + ${MPIEXEC_PREFLAGS} + ./${testname} + ${MPIEXEC_POSTFLAGS} + ) + elseif(${LIBNEGF_TEST_WITH_MPIEXEC} AND NOT ${LIBNEGF_TEST_MPI}) + set(num_procs 1) + set(cmd ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${num_procs} ${MPIEXEC_PREFLAGS} + ./${testname} ${MPIEXEC_POSTFLAGS} + ) + else() + set(cmd ${testname}) + endif() + add_test(NAME ${testname} COMMAND ${cmd}) + set_property(TEST ${testname} PROPERTY LABELS ${LIBNEGF_TEST_LABELS}) + set_property(TEST ${testname} PROPERTY PROCESSORS ${num_procs}) + + transfer_test_data(${testname}) +endfunction() foreach(test-directory IN LISTS test-directories) - add_subdirectory(${test-directory}) + add_subdirectory(${test-directory}) endforeach() + +add_executable(ext_system_test ext_system_test.f90) +target_link_libraries(ext_system_test syscalls_objlib) +add_test(ext_system_test ext_system_test) diff --git a/tests/c_int/CMakeLists.txt b/tests/c_int/CMakeLists.txt index 944afcef..56f40f3a 100644 --- a/tests/c_int/CMakeLists.txt +++ b/tests/c_int/CMakeLists.txt @@ -1 +1 @@ -setup_c_test(c_int "normal" hello.c) +setup_test(c_int hello.c) diff --git a/tests/c_int_elph_deph/CMakeLists.txt b/tests/c_int_elph_deph/CMakeLists.txt index c277efad..4a25b6a5 100644 --- a/tests/c_int_elph_deph/CMakeLists.txt +++ b/tests/c_int_elph_deph/CMakeLists.txt @@ -1 +1 @@ -setup_c_test(c_int_elph_deph "normal" c_int_elph_deph.c) +setup_test(c_int_elph_deph c_int_elph_deph.c) diff --git a/tests/c_int_file/CMakeLists.txt b/tests/c_int_file/CMakeLists.txt index 97d0d72c..bd2b0526 100644 --- a/tests/c_int_file/CMakeLists.txt +++ b/tests/c_int_file/CMakeLists.txt @@ -1 +1 @@ -setup_c_test(c_int_file "normal" hello.c) +setup_test(c_int_file hello.c) diff --git a/tests/c_mpi_init/CMakeLists.txt b/tests/c_mpi_init/CMakeLists.txt index a81a8ddb..c643a2c0 100644 --- a/tests/c_mpi_init/CMakeLists.txt +++ b/tests/c_mpi_init/CMakeLists.txt @@ -1 +1 @@ -setup_c_test(c_mpi_init "normal" c_mpi_init.c) +setup_test(c_mpi_init c_mpi_init.c) diff --git a/tests/cpp_mpi_init/CMakeLists.txt b/tests/cpp_mpi_init/CMakeLists.txt index e349010d..16250a9a 100644 --- a/tests/cpp_mpi_init/CMakeLists.txt +++ b/tests/cpp_mpi_init/CMakeLists.txt @@ -1 +1 @@ -setup_cpp_test(cpp_mpi_init "normal" cpp_mpi_init.cpp) +setup_test(cpp_mpi_init cpp_mpi_init.cpp MPI) diff --git a/tests/cpp_mpi_init/cpp_mpi_init.cpp b/tests/cpp_mpi_init/cpp_mpi_init.cpp index 6ddc4ef0..9d3c1835 100644 --- a/tests/cpp_mpi_init/cpp_mpi_init.cpp +++ b/tests/cpp_mpi_init/cpp_mpi_init.cpp @@ -1,37 +1,52 @@ - +/*!!--------------------------------------------------------------------------! + *!! libNEGF: a general library for Non-Equilibrium Greens functions. ! + *!! Copyright (C) 2012 - 2026 ! + *!! ! + *!! This file is part of libNEGF: a library for ! + *!! Non Equilibrium Green's Function calculation ! + *!! ! + *!! Developers: Alessandro Pecchia, Daniele Soccodato ! + *!! Former Contributors: Gabriele Penazzi, Luca Latessa, Aldo Di Carlo ! + *!! ! + *!! libNEGF is free software: you can redistribute it and/or modify ! + *!! it under the terms of the GNU Lesse General Public License as published ! + *!! by the Free Software Foundation, either version 3 of the License, or ! + *!! (at your option) any later version. ! + *!! ! + *!! You should have received a copy of the GNU Lesser General Public ! + *!! License along with libNEGF. If not, see ! + *!! . ! + *!!--------------------------------------------------------------------------! + */ #include "libnegf.hpp" -#include -#include + #include -int main() -{ - int handler[NEGF_HSIZE]; - int *hand = &handler[0]; - int ierr; +#include +#include - ierr = MPI_Init(NULL, NULL); - if (ierr != 0) - { - printf("Error in mpi_init \n"); - return 1; - } - printf("Initializing libNEGF \n"); - negf_init_session(hand); - negf_init(hand); +int main() { + int ierr = MPI_Init(nullptr, nullptr); + if(ierr != 0) { + fprintf(stderr, "Error in mpi_init: %d\n", ierr); + return 1; + } - MPI_Fint global_comm_f = MPI_Comm_c2f(MPI_COMM_WORLD); - negf_set_mpi_fcomm(hand, global_comm_f); - MPI_Fint cart_comm, k_comm, en_comm; - negf_cartesian_init(hand, global_comm_f, 1, cart_comm, k_comm, en_comm); + printf("Initializing libNEGF\n"); + int handler[NEGF_HSIZE]; + negf_init_session(handler); + negf_init(handler); - //Release library - negf_destruct_libnegf(hand); - negf_destruct_session(hand); - printf("Done \n"); + MPI_Fint global_comm_f = MPI_Comm_c2f(MPI_COMM_WORLD); + negf_set_mpi_fcomm(handler, global_comm_f); + MPI_Fint cart_comm, k_comm, en_comm; + negf_cartesian_init(handler, global_comm_f, 1, cart_comm, k_comm, en_comm); - MPI_Finalize(); + //Release library + negf_destruct_libnegf(handler); + negf_destruct_session(handler); + printf("Done\n"); - return 0; + MPI_Finalize(); } diff --git a/ext_system/test.f90 b/tests/ext_system_test.f90 similarity index 66% rename from ext_system/test.f90 rename to tests/ext_system_test.f90 index 1f30d2e6..034c104c 100644 --- a/ext_system/test.f90 +++ b/tests/ext_system_test.f90 @@ -4,7 +4,7 @@ program test character(100) :: folder - folder = "testfolder" + folder = "ext_system_testfolder" call create_directory(trim(folder)) @@ -12,9 +12,9 @@ program test write(101, *) 'test', 101 close(101) - !call remove_file(trim(folder)//'/afile') + call remove_file(trim(folder)//'/afile') - !call remove_directory(folder) + call remove_directory(folder) end program test diff --git a/tests/f90Si2x2/CMakeLists.txt b/tests/f90Si2x2/CMakeLists.txt index 840a2cec..57f108db 100644 --- a/tests/f90Si2x2/CMakeLists.txt +++ b/tests/f90Si2x2/CMakeLists.txt @@ -1,10 +1,7 @@ -set(sources - constants.F90 - readHS.F90 - matconv.F90 - test.F90) - -setup_f90_mpi_test(f90Si2x2 "normal" ${sources}) - +setup_test( + f90Si2x2 test.F90 + MPI + SOURCES constants.F90 matconv.F90 readHS.F90 + LABELS long +) unzip_test_data(f90Si2x2 "hs.tar.xz") - diff --git a/tests/f90Si_nin/CMakeLists.txt b/tests/f90Si_nin/CMakeLists.txt index 7721c875..69b6034c 100644 --- a/tests/f90Si_nin/CMakeLists.txt +++ b/tests/f90Si_nin/CMakeLists.txt @@ -1,10 +1,7 @@ -set(sources - constants.F90 - readHS.F90 - matconv.F90 - test.F90) - -setup_f90_mpi_test(f90Si_nin "long" ${sources}) - +setup_test( + f90Si_nin test.F90 + MPI + SOURCES constants.F90 matconv.F90 readHS.F90 + LABELS long +) unzip_test_data(f90Si_nin "hs.tar.xz") - diff --git a/tests/f90Si_nin/transmission_ref.dat b/tests/f90Si_nin/transmission_ref.dat index afc74509..39516429 100644 --- a/tests/f90Si_nin/transmission_ref.dat +++ b/tests/f90Si_nin/transmission_ref.dat @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:020e400bdde24a9fdf260ea757b134afda07bda65b725b62cfb3a86bad60ad30 +oid sha256:20d144dba082e78dae2d676ae8871812d4ad7aa37879a44ff2de9f8d7424677f size 784 diff --git a/tests/f90Si_nin_40pl/CMakeLists.txt b/tests/f90Si_nin_40pl/CMakeLists.txt index a1643d9a..5a79e652 100644 --- a/tests/f90Si_nin_40pl/CMakeLists.txt +++ b/tests/f90Si_nin_40pl/CMakeLists.txt @@ -1,10 +1,7 @@ -set(sources - constants.F90 - readHS.F90 - matconv.F90 - test.F90) - -setup_f90_mpi_test(f90Si_nin_40pl "long" ${sources}) - +setup_test( + f90Si_nin_40pl test.F90 + MPI + SOURCES constants.F90 matconv.F90 readHS.F90 + LABELS long +) unzip_test_data(f90Si_nin_40pl "hs.tar.xz") - diff --git a/tests/f90Si_nin_40pl/test.F90 b/tests/f90Si_nin_40pl/test.F90 index 78b2c7fa..ba8f3839 100644 --- a/tests/f90Si_nin_40pl/test.F90 +++ b/tests/f90Si_nin_40pl/test.F90 @@ -30,8 +30,16 @@ program test call mpifx_init(impierr); call globalComm%init(); - nGroups = 1 + if (globalComm%size < 4) then + if (globalComm%lead) then + write(*,*) 'ERROR: this test is designed to run on at least 4 processes' + end if + stop + else + nGroups = 4 + end if + pnegf => negf if (globalComm%lead) write(*,*) 'Initializing libNEGF' @@ -127,10 +135,20 @@ program test ! Here we set the parameters, only the ones different from default ! ---------------------------------------------------------------------------------- call get_params(pnegf, params) - params%verbose = 100 - params%Emin = mu(1)-10.0_dp*kt(1)+0.012_dp ! -0.298 - params%Emax = mu(2)+10.0_dp*kt(1)-0.010_dp ! +0.30 - params%Estep = 0.002_dp + params%verbose = 100 + ! The following parameters are good for benchmarks (300 Epoints) + !params%Emin = mu(1)-10.0_dp*kt(1)+0.012_dp/eV2Hartree ! -0.298 + !params%Emax = mu(2)+10.0_dp*kt(1)-0.010_dp/eV2Hartree ! +0.30 + !params%Estep = 0.002_dp/eV2Hartree + + ! The following parameters are used for testing (20 Epoints) + params%Emax = mu(1)+10.0_dp*kt(1) ! mu+0.26 + params%Emin = mu(1)-0.12_dp/eV2Hartree ! mu-0.13 + params%Estep = 0.02_dp/eV2Hartree + + if (globalComm%lead) then + write(*,*) 'nPoints=',(params%Emax-params%Emin)/params%Estep + 1 + end if ! nStep = 0.598/0.002 + 1 = 300 params%delta = 1.e-5_dp !Already in Hartree params%mu(1:2) = mu diff --git a/tests/f90Si_nin_40pl/transmission_ref.dat b/tests/f90Si_nin_40pl/transmission_ref.dat new file mode 100644 index 00000000..ca415a19 --- /dev/null +++ b/tests/f90Si_nin_40pl/transmission_ref.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:960ada2bac077c2f149093643bd7d6c1b2a5c1a12580f019817e6bab676155c1 +size 980 diff --git a/tests/f90_transmission/CMakeLists.txt b/tests/f90_transmission/CMakeLists.txt index ff77843b..379acd57 100644 --- a/tests/f90_transmission/CMakeLists.txt +++ b/tests/f90_transmission/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(f90_transmission "normal" hello.F90) +setup_test(f90_transmission hello.F90) diff --git a/tests/f90elph_deph/CMakeLists.txt b/tests/f90elph_deph/CMakeLists.txt index 3be2bbea..27a405da 100644 --- a/tests/f90elph_deph/CMakeLists.txt +++ b/tests/f90elph_deph/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(f90elph_deph "normal" hello.F90) +setup_test(f90elph_deph hello.F90) diff --git a/tests/f90init/CMakeLists.txt b/tests/f90init/CMakeLists.txt index e1421948..928750fb 100644 --- a/tests/f90init/CMakeLists.txt +++ b/tests/f90init/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(f90init "normal" f90init.F90) +setup_test(f90init f90init.F90) diff --git a/tests/f90int/CMakeLists.txt b/tests/f90int/CMakeLists.txt index 40128108..4b9b66ff 100644 --- a/tests/f90int/CMakeLists.txt +++ b/tests/f90int/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(f90int "normal" hello.F90) +setup_test(f90int hello.F90) diff --git a/tests/f90int_file/CMakeLists.txt b/tests/f90int_file/CMakeLists.txt index 714a6da7..c10406b4 100644 --- a/tests/f90int_file/CMakeLists.txt +++ b/tests/f90int_file/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(f90int_file "normal" hello.F90) +setup_test(f90int_file hello.F90) diff --git a/tests/f90read_hs/CMakeLists.txt b/tests/f90read_hs/CMakeLists.txt index 888f10a4..3d37ea20 100644 --- a/tests/f90read_hs/CMakeLists.txt +++ b/tests/f90read_hs/CMakeLists.txt @@ -1,7 +1,5 @@ -set(sources - constants.F90 - readHS.F90 - matconv.F90 - test.F90) - -setup_f90_test(f90read_hs "normal unit" ${sources}) +setup_test( + f90read_hs test.F90 + SOURCES constants.F90 matconv.F90 readHS.F90 + LABELS unit +) diff --git a/tests/f90read_hs/readHS.F90 b/tests/f90read_hs/readHS.F90 index 7c3a10b6..28ed4adb 100644 --- a/tests/f90read_hs/readHS.F90 +++ b/tests/f90read_hs/readHS.F90 @@ -415,7 +415,7 @@ subroutine apply_shifts(ham, over, shift) real(dp) :: tmpH(orb%mOrb,orb%mOrb), tmpS(orb%mOrb,orb%mOrb) real(dp) :: shiftBlk1(orb%mOrb,orb%mOrb) real(dp) :: shiftBlk2(orb%mOrb,orb%mOrb) - + nAtom = size(orb%nOrbAtom) shiftBlk1 = 0.0_dp shiftBlk2 = 0.0_dp @@ -427,8 +427,8 @@ subroutine apply_shifts(ham, over, shift) iAt2f = img2CentCell(iAt2) nOrb2 = orb%nOrbAtom(iAt2f) iOrig1 = iPair(iNeigh, iAt1) - iBlk1 = iAtomStart(iAt1)-1 - iBlk2 = iAtomStart(iAt2f)-1 + iBlk1 = iAtomStart(iAt1-1)-1 + iBlk2 = iAtomStart(iAt2f-1)-1 do ii = 1, nOrb1 shiftBlk1(ii,ii) = shift(iBlk1+ii) end do diff --git a/tests/testCUDA/CMakeLists.txt b/tests/testCUDA/CMakeLists.txt index 5c8bfcef..08114436 100644 --- a/tests/testCUDA/CMakeLists.txt +++ b/tests/testCUDA/CMakeLists.txt @@ -1 +1 @@ -setup_f90_test(testCUDA "normal cuda" test.f90) +setup_test(testCUDA test.f90 GPUONLY) diff --git a/tests/testCUDA/test.f90 b/tests/testCUDA/test.f90 index 069624c2..fbff1871 100644 --- a/tests/testCUDA/test.f90 +++ b/tests/testCUDA/test.f90 @@ -1,7 +1,8 @@ program test use ln_precision + use ln_allocation use mat_def - use lib_param, only : cublasHandle + use lib_param, only : cublasHandle, cusolverDnHandle use libnegf, only : cublasInitialize, cublasFinalize use cudautils implicit none @@ -11,29 +12,48 @@ program test !end type TArray integer, parameter :: Npl = 5 - integer, parameter :: Nr = 2000, Nc = 3000 + integer, parameter :: Nr = 5000, Nc = 300 type(z_DNS), allocatable :: A(:,:) - type(z_DNS) :: work1, work2, myMat, myMat_dag + type(z_DNS) :: work1, work2, work3 complex(dp) :: alpha, beta real(dp), allocatable :: R(:,:) - integer :: istat - type(cublasHandle) :: hcublas + real(dp) :: tr1, tr2 + integer :: istat, ii, ndevs + integer(8) :: totalmem, freemem + type(cublasHandle) :: hblas + type(cusolverDnHandle) :: hsolver - call cublasInitialize(hcublas) + call getDeviceCount(ndevs) + print*,'found ',ndevs,' GPUs' + call getDevice(ii) + print*,'current device number: ',ii + + call getDevMemInfo(freemem, totalmem) + print*, 'Total device memory: ',totalmem + print*, ' Free device memory: ',freemem + + call cublasInitialize(hblas) + call cusolverInitialize(hsolver) allocate(A(Npl,Npl)) print*,'create A11' call createAll(A(1,1), Nr, Nr) print*,'create A12' call createAll(A(1,2), Nr, Nc) + print*,'create A22' + call createAll(A(2,2), Nr, Nr) allocate(R(Nr,Nr)) call random_number(R) A(1,1)%val = R + A(2,2)%val = R call random_number(R) A(1,1)%val = A(1,1)%val + (0.0_dp,1.0_dp)*R + A(2,2)%val = A(2,2)%val + (0.0_dp,1.0_dp)*R deallocate(R) - + print*,'maxval abs(A11)',maxval(abs(A(1,1)%val)) + + allocate(R(Nr,Nc)) call random_number(R) A(1,2)%val = R @@ -41,91 +61,288 @@ program test A(1,2)%val = A(1,2)%val + (0.0_dp,1.0_dp)*R deallocate(R) - !print*,'createGPU A11' - !call createGPU(A(1,1)) - print*,'copyGPU A11' + !///////////////////////////////////////////////////////////////////// + istat = 0 + + print*, '' + print*,'1.1 TEST matmul GPU: C = A*B' + print*,'copy A11 to GPU' call copyToGPU(A(1,1)) - !print*,'createGPU A12' - !call createGPU(A(1,2)) print*,'copyGPU A12' call copyToGPU(A(1,2)) print*,'create work' call createAll(work1, Nr, Nc) - !print*,'createGPU work' - !call createGPU(work1) - print*,'matmul' alpha = cmplx(1.0,0.0,dp) beta = cmplx(0.0,0.0,dp) - call matmul_gpu(hcublas,alpha,A(1,1),A(1,2),beta,work1) + call matmul_gpu(hblas,alpha,A(1,1),A(1,2),beta,work1) print*,'copy Back' call copyFromGPU(work1) + print*,'matmul CPU' call create(work2, Nr, Nc) - work2%val = matmul(A(1,1)%val, A(1,2)%val) + print*,'TEST check' if (any(abs(work1%val-work2%val) > 1e-10)) then print*,'TEST failed' + istat = 1001 else print*,'TEST Ok' end if + print*,'-------------------------' - !print*, work%val(1,1), work%val(2,1), work%val(3,1) - print*,'finalize' - !call deleteGPU(work1) - !call deleteGPU(A(1,2)) - - call destroyAll(work1) + call createAll(A(2,1), Nc, Nr) + A(2,1)%val = conjg(transpose(A(1,2)%val)) call destroyAll(A(1,2)) + call copyToGPU(A(2,1)) + + print*,'1.1 TEST matmul GPU: C = A*B^H' + call matmul_gpu(hblas,alpha,A(1,1),A(2,1),beta,work1,'dag_2nd') + + print*,'copy Back' + call copyFromGPU(work1) + + print*,'TEST check' + if (any(abs(work1%val-work2%val) > 1e-10)) then + print*,'TEST failed' + istat = 10021 + else + print*,'TEST Ok' + end if + + print*,'1.2 TEST matmul GPU: C = C + A*B^H' + beta = cmplx(1.0,0.0,dp) + call matmul_gpu(hblas,alpha,A(1,1),A(2,1),beta,work1,'dag_2nd') + + print*,'copy Back' + call copyFromGPU(work1) + + print*,'TEST check' + if (any(abs(work1%val-2.0*work2%val) > 1e-10)) then + print*,'TEST failed' + istat = 10022 + else + print*,'TEST Ok' + end if + + call destroyAll(A(2,1)) + call destroyAll(work1) call destroy(work2) + print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// - print*, '' - print*, 'TEST dagger' + print*,'3. TEST trace' + call trace_gpu(hblas, A(1,1), tr1) + print*,'TEST check' + tr2 = real(0.0,dp) + do ii = 1, Nr + tr2 = tr2 + real(A(1,1)%val(ii,ii)) + end do + if (abs(tr1-tr2) > 1e-10) then + print*,'TEST failed', tr1, tr2 + istat = 1003 + else + print*,'TEST Ok' + end if + print*,'-------------------------' + + print* + print*,'4. TEST hardcopy mat ' + print*,'copy A22 to GPU' + call copyToGPU(A(2,2)) + + print*,'create work' call createAll(work1, Nr, Nr) - call dagger_gpu(hcublas, A(1,1), work1) + + !///////////////////////////////////////////////////////////////////// + call copy_mat_gpu(hblas, A(2,2), work1) + print*,'copy Back' + call copyFromGPU(work1) + print*,'TEST check' + if (any(abs(work1%val-A(2,2)%val) > 1e-13)) then + print*,'TEST failed' + istat = 1004 + else + print*,'TEST Ok' + end if + print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// + print*, '' + print*, '5. TEST dagger GPU' + call dagger_gpu(hblas, A(1,1), work1) + print*,'copy Back' call copyFromGPU(work1) + print*, 'dagger CPU' call create(work2, Nr, Nr) work2%val = conjg(transpose(A(1,1)%val)) + print*,'TEST check' if (any(abs(work1%val-work2%val) > 1e-10)) then print*,'TEST failed' + istat = 1005 + else + print*,'TEST Ok' + end if + print*,'finalize' + call destroy(work2) + print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// + !print* + !print*,'6. TEST kernelsum GPU' + !alpha = cmplx(1.0,0.0,dp) + !beta = cmplx(1.0,0.0,dp) + !call kernelsum_gpu(work1, alpha, A(1,1), beta, A(2,2)) + !print*,'copy Back' + !call copyFromGPU(work1) + !print*,'TEST check' + !if (any(abs(work1%val-A(1,1)%val-A(2,2)%val) > 1e-10)) then + ! print*,'TEST failed' + ! istat = 1006 + !else + ! print*,'TEST Ok' + !end if + !print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// + print* + print*,'7. TEST matsum GPU C = A + B' + alpha = cmplx(1.0,0.0,dp) + beta = cmplx(1.0,0.0,dp) + call matsum_gpu(hblas, alpha, A(1,1), beta, A(2,2), work1) + print*,'copy Back' + call copyFromGPU(work1) + + print*,'TEST check' + if (any(abs(work1%val-A(1,1)%val-A(2,2)%val) > 1e-10)) then + print*,'TEST failed' + istat = 1007 + else + print*,'TEST Ok' + end if + print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// + print* + print*,'7b. TEST matsum in-place GPU C = C + B' + call matsum_gpu(hblas, alpha, work1, beta, A(2,2), work1) + print*,'copy Back' + call copyFromGPU(work1) + print*,'TEST check' + if (any(abs(work1%val-A(1,1)%val-2.0_dp*A(2,2)%val) > 1e-10)) then + print*,'TEST failed' + istat = 1007 + else + print*,'TEST Ok' + end if + call destroyAll(A(2,2)) + print*,'-------------------------' + + + !///////////////////////////////////////////////////////////////////// + print*, '8. TEST spectral as C = i[A - A^H]' + alpha = cmplx(0.0,1.0,dp) + beta = cmplx(0.0,-1.0,dp) + call matsum_gpu(hblas, alpha, A(1,1), beta, A(1,1), work1,'dag_2nd') + print*,'copy Back' + call copyFromGPU(work1) + print*,'TEST check' + if (any(abs(A(1,1)%val-conjg(transpose(A(1,1)%val))+alpha*work1%val) > 1e-10)) then + print*,'TEST failed' + istat = 1008 + else + print*,'TEST Ok' + end if + print*,'finalize' + print*,'-------------------------' + + + print* + print*,'9. TEST spectral GPU' + call spectral_gpu(hblas, A(1,1), work1) + + print*,'copy Back' + call copyFromGPU(work1) + + print*,'TEST check' + beta = cmplx(0.0,1.0,dp) + if (any(abs(A(1,1)%val-conjg(transpose(A(1,1)%val))+beta*work1%val) > 1e-10)) then + print*,'TEST failed' + istat = 1009 else print*,'TEST Ok' end if print*,'finalize' - call destroyAll(work1) call destroyAll(A(1,1)) - call destroy(work2) + print*,'-------------------------' + + !///////////////////////////////////////////////////////////////////// + print*, '' + print*, '10. TEST inverse GPU' + print*,'create A11' + call createAll(A(1,1), Nr, Nr) + allocate(R(Nr,Nr)) + call random_number(R) + A(1,1)%val = R + call random_number(R) + A(1,1)%val = A(1,1)%val + (0.0_dp,1.0_dp)*R + deallocate(R) + print*,'maxval abs(A11)',maxval(abs(A(1,1)%val)) - !call create(myMat, 4, 4) - !myMat%val(:,:) = (0.0_dp, 0.0_dp) - !myMat%val(1,1) = (1.0_dp,1.0_dp) - !myMat%val(1,2) = (0.0_dp,-1.0_dp) - !myMat%val(1,3) = (0.0_dp,4.0_dp) - !myMat%val(1,4) = (3.0_dp,0.0_dp) - !print*, 'A:' - !call print_myMat(myMat%val) - !print*, '' - - !call copyToGPU(myMat) - !call create(myMat_dag, 4, 4) - !call createGPU(myMat_dag) - - !call copyFromGPU(myMat_dag) - !call print_myMat(myMat_dag%val) + call copyToGPU(A(1,1)) + call createAll(work1, Nr, Nr) + + call inverse_gpu(hblas, hsolver, A(1,1), work1, istat) + if (istat > 0) then + print*, 'Error in inverse:',istat + istat = 1010 + end if + + print*,'TEST A*A^-1 = Id' + call createAll(work2, Nr, Nr) + alpha = cmplx(1.0,0.0,dp) + beta = cmplx(0.0,0.0,dp) + call matmul_gpu(hblas,alpha,A(1,1),work1,beta,work2) + print*,'copy Back' + call copyFromGPU(work2) - !call deleteGPU(myMat) - !call deleteGPU(myMat_dag) - !call destroy(myMat) - !call destroy(myMat_dag) + ! Set an identity matrix + work1%val = (0.0_dp, 0.0_dp) + do ii = 1, Nr + work1%val(ii,ii) = (1.0_dp, 0.0_dp) + end do + + print*,'check' + if (any(abs(work1%val-work2%val) > 1e-10)) then + print*,'TEST failed' + istat = 1010 + else + print*,'TEST Ok' + end if - call cublasFinalize(hcublas) + print*,'finalize' + call destroyAll(work1) + call destroyAll(A(1,1)) + call destroyAll(work2) + print*,'-------------------------' + print*,"Memory summary:" + call writePeakInfo(6) + call writeMemInfo(6) + + call cublasFinalize(hblas) + call cusolverFinalize(hsolver) + + if (istat /= 0) then + error stop "Error found in one of the HIP tests" + end if contains subroutine print_myMat(mat) diff --git a/tests/testCUDA_decimation/CMakeLists.txt b/tests/testCUDA_decimation/CMakeLists.txt index 5e542596..5d1c4442 100644 --- a/tests/testCUDA_decimation/CMakeLists.txt +++ b/tests/testCUDA_decimation/CMakeLists.txt @@ -1,5 +1,5 @@ -set(sources - main.f90 - random.f90) - -setup_f90_test(testCUDA_decimation "normal cuda" ${sources}) +setup_test( + testCUDA_decimation main.f90 + GPUONLY + SOURCES random.f90 +) diff --git a/tests/testCUDA_decimation/main.f90 b/tests/testCUDA_decimation/main.f90 index bea41544..73b97bc1 100644 --- a/tests/testCUDA_decimation/main.f90 +++ b/tests/testCUDA_decimation/main.f90 @@ -2,14 +2,14 @@ program test1 use ln_precision use libnegf use random - use contselfenergy, only : decimation_gpu + use cudautils, only : decimation_gpu implicit none complex(dp), dimension(:,:), allocatable :: H00, H01 complex(dp), dimension(:,:), allocatable :: S00, S01 complex(dp), dimension(:,:), allocatable, target :: A0, B0, C0, G0 character(10) :: arg integer :: N, ncyc, fu, t1, t2, cr, cm, stat - real(dp) :: delta, En + real(dp) :: delta, En, acc=1e-12_dp complex(dp) :: z logical :: readfiles @@ -75,10 +75,10 @@ program test1 call init_negf(negf) call system_clock(t1, cr, cm) - call decimation_gpu(negf, G0, A0, B0, C0, N, .false., ncyc) + call decimation_gpu(negf, G0, A0, B0, C0, N, ncyc, acc, loginfo=.true.) call system_clock(t2, cr, cm) - write(*,*) "decimation_cpu time: ",(t2-t1)*1.0/cr,"sec" + write(*,*) "decimation_gpu time: ",(t2-t1)*1.0/cr,"sec" print*,'converged in',ncyc,' iterations' deallocate(A0,B0,C0,G0)