Draft release 0.2

.github/workflows/cmake-cuda-multi-compiler.yml

@@ -21,7 +21,7 @@ jobs:
         cpp_compiler: [g++, clang++]
 
     steps:
-      - uses: Jimver/cuda-toolkit@v0.2.15
+      - uses: Jimver/cuda-toolkit@v0.2.22
         id: cuda-toolkit
         with:
           method: 'network'
@@ -46,6 +46,9 @@ jobs:
           -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
           -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON
           -DSPDLOG_FMT_EXTERNAL:BOOL=ON
+          -DCOMPAS_USE_CUDA:BOOL=ON
+          -DCOMPAS_BUILD_TESTS:BOOL=ON
+          -DCOMPAS_BUILD_BENCHMARKS:BOOL=ON
           -S ${{ github.workspace }}
 
       - name: Build

.github/workflows/cmake-hip.yml

@@ -0,0 +1,58 @@
+# This starter workflow is for a CMake project running on multiple platforms. There is a different starter workflow if you just want a single platform.
+# See: https://github.com/actions/starter-workflows/blob/main/ci/cmake-single-platform.yml
+name: HIP build
+
+on:
+  push:
+    branches: [ "dev" ]
+  pull_request:
+    branches: [ "main" ]
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+
+    strategy:
+      fail-fast: true
+
+      matrix:
+        os: [ubuntu-latest]
+        build_type: [Debug]
+
+    steps:
+      - uses: loostrum/rocm-installer@main
+        with:
+          version: 6.2.2
+          packages: hip-dev hipcc rocm-device-libs rocblas rocm-hip-runtime-dev rocrand
+
+      - uses: actions/checkout@v4
+        with:
+          submodules: 'recursive'
+
+      - name: Set reusable strings
+        # Turn repeated input strings (such as the build output directory) into step outputs. These step outputs can be used throughout the workflow file.
+        id: strings
+        shell: bash
+        run: |
+          echo "build-output-dir=${{ github.workspace }}/build" >> "$GITHUB_OUTPUT"
+
+      - name: Configure CMake
+        # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+        # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        run: >
+          cmake -B ${{ steps.strings.outputs.build-output-dir }}
+          -DCMAKE_PREFIX_PATH=/opt/rocm
+          -DCMAKE_HIP_COMPILER_ROCM_ROOT:PATH=/opt/rocm
+          -DCMAKE_HIP_ARCHITECTURES=gfx90a
+          -DAMDGPU_TARGETS=gfx90a
+          -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
+          -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON
+          -DSPDLOG_FMT_EXTERNAL:BOOL=ON
+          -DCOMPAS_USE_HIP:BOOL=ON
+          -DCOMPAS_BUILD_TESTS:BOOL=ON
+          -DCOMPAS_BUILD_BENCHMARKS:BOOL=ON
+          -S ${{ github.workspace }}
+
+      - name: Build
+        # Build your program with the given configuration. Note that --config is needed because the default Windows generator is a multi-config generator (Visual Studio generator).
+        run: cmake --build ${{ steps.strings.outputs.build-output-dir }} --config ${{ matrix.build_type }} --parallel 4

.gitmodules

@@ -4,6 +4,9 @@
 [submodule "thirdparty/kernel_launcher"]
 	path = thirdparty/kernel_launcher
 	url = https://github.com/KernelTuner/kernel_launcher.git
+[submodule "thirdparty/kernel_float"]
+	path = thirdparty/kernel_float
+    url = https://github.com/KernelTuner/kernel_float.git
 [submodule "thirdparty/kmm"]
 	path = thirdparty/kmm
-	url = https://github.com/NLeSC-COMPAS/kmm.git
+	url = http://github.com/NLeSC-COMPAS/kmm.git

CMakeLists.txt

@@ -1,7 +1,31 @@
 cmake_minimum_required(VERSION 3.10)
+set(CMAKE_POLICY_DEFAULT_CMP0077 NEW)
 
 set(PROJECT_NAME "compas-toolkit")
-project(${PROJECT_NAME} LANGUAGES CXX CUDA)
+project(${PROJECT_NAME} LANGUAGES CXX VERSION 0.2)
+
+# User options (features)
+option(COMPAS_USE_CUDA "Use the CUDA backend (default)" OFF)
+option(COMPAS_USE_HIP  "Use the HIP backend" OFF)
+option(COMPAS_BUILD_TEST "Build tests" OFF)
+option(COMPAS_BUILD_BENCHMARK "Build benchmarks" OFF)
+option(COMPAS_ENABLE_LINTER  "Enable clang-tidy linter" OFF)
+
+# If neither backend is enabled, default to CUDA
+if(NOT COMPAS_USE_CUDA AND NOT COMPAS_USE_HIP)
+    set(COMPAS_USE_CUDA ON CACHE BOOL "Use the CUDA backend (default)" FORCE)
+    message(STATUS "No GPU backend was enabled. Compiling for CUDA.")
+endif()
+
+# Check for mutual exclusivity
+if(COMPAS_USE_CUDA AND COMPAS_USE_HIP)
+    message(FATAL_ERROR "Exactly one backend between CUDA and HIP must be enabled.")
+endif()
+
+# Set CMAKE_CUDA_ARCHITECTURES if not defined. This needs to happen before calling `add_library`
+if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
+    set(CMAKE_CUDA_ARCHITECTURES all-major)
+endif()
 
 # Enable C++17 support
 set(CXX_STANDARD 17)
@@ -21,29 +45,71 @@ install(TARGETS ${PROJECT_NAME}
 target_include_directories(${PROJECT_NAME} PUBLIC "${PROJECT_SOURCE_DIR}/include")
 target_include_directories(${PROJECT_NAME} PUBLIC "${PROJECT_SOURCE_DIR}/src")
 
-# The CUDA compiler struggles with the "Werror" options, so we need to explicitly forward it to the host compiler
-# using `-Xcompiler=-Werror` if we are compiling CUDA code.
-set(CXXFLAGS
-        $<$<COMPILE_LANGUAGE:CUDA>:-forward-unknown-to-host-compiler>
-        $<$<COMPILE_LANGUAGE:CUDA>:--generate-line-info>
-        -Wall -Wextra -Wconversion -Wno-unused-parameter
-        $<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=-Werror>
-        $<$<COMPILE_LANGUAGE:CUDA>:--use_fast_math>
-        $<$<COMPILE_LANGUAGE:CUDA>:-Xptxas="-v">
-)
-target_compile_options(${PROJECT_NAME} PRIVATE ${CXXFLAGS})
-
 # Enable PIC
 set_property(TARGET ${PROJECT_NAME} PROPERTY POSITION_INDEPENDENT_CODE ON)
 
-set(PROJECT_CLANG_TIDY clang-tidy)
-set_target_properties(${PROJECT_NAME} PROPERTIES CXX_CLANG_TIDY "${PROJECT_CLANG_TIDY}")
+if(COMPAS_ENABLE_LINTER)
+    set(PROJECT_CLANG_TIDY clang-tidy)
+    set_target_properties(${PROJECT_NAME} PROPERTIES CXX_CLANG_TIDY "${PROJECT_CLANG_TIDY}")
+endif()
+
+if(COMPAS_USE_CUDA)
+    enable_language(CUDA)
+    set(CUDA_SEPARABLE_COMPILATION ON)
+    set(CMAKE_CUDA_SEPARABLE_COMPILATION ON)
+    set(CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+    set_property(TARGET ${PROJECT_NAME} PROPERTY CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+
+    find_package(CUDAToolkit REQUIRED)
+    target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cudart_static)
+    target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cuda_driver)
+    target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cublas)
+    target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::nvrtc)
+    set_target_properties(${PROJECT_NAME} PROPERTIES CUDA_ARCHITECTURES ${CMAKE_CUDA_ARCHITECTURES})
+
+    # The CUDA compiler struggles with the "Werror" options, so we need to explicitly forward it to the host compiler
+    # using `-Xcompiler=-Werror` if we are compiling CUDA code.
+    set(CXXFLAGS
+            $<$<COMPILE_LANGUAGE:CUDA>:-forward-unknown-to-host-compiler>
+            $<$<COMPILE_LANGUAGE:CUDA>:--generate-line-info>
+            -Wall -Wextra -Wconversion -Wno-unused-parameter
+            #$<$<COMPILE_LANGUAGE:CUDA>:-Xcompiler=-Werror>
+            $<$<COMPILE_LANGUAGE:CUDA>:--use_fast_math>
+            $<$<COMPILE_LANGUAGE:CUDA>:-Xptxas="-v">
+    )
+    target_compile_options(${PROJECT_NAME} PRIVATE ${CXXFLAGS})
 
-find_package(CUDAToolkit REQUIRED)
-target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cudart_static)
-target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cuda_driver)
-target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::cublas)
-target_link_libraries(${PROJECT_NAME} PUBLIC CUDA::nvrtc)
+    # Define `COMPAS_USE_CUDA` macro so that headers can detect CUDA usage
+    target_compile_definitions(${PROJECT_NAME} PUBLIC COMPAS_USE_CUDA=1)
+    set(KMM_USE_CUDA ON)
+elseif(COMPAS_USE_HIP)
+    if(NOT DEFINED HIP_PATH)
+        if(NOT DEFINED ENV{HIP_PATH})
+            set(HIP_PATH "/opt/rocm/hip" CACHE PATH "Path to which HIP has been installed")
+        else()
+            set(HIP_PATH $ENV{HIP_PATH} CACHE PATH "Path to which HIP has been installed")
+        endif()
+    endif()
+    set(CMAKE_MODULE_PATH "${HIP_PATH}/cmake" ${CMAKE_MODULE_PATH})
+    enable_language(HIP)
+    set_source_files_properties(${sources} PROPERTIES LANGUAGE HIP)
+
+    find_package(HIP REQUIRED)
+    find_package(ROCBLAS REQUIRED)
+    target_link_libraries(${PROJECT_NAME} PUBLIC hip::host)
+    target_link_libraries(${PROJECT_NAME} PUBLIC hip::device)
+    target_link_libraries(${PROJECT_NAME} PUBLIC roc::rocblas)
+
+    # Enable HIP in KMM and kernel float
+    set(KMM_USE_HIP ON)
+    set(KERNEL_FLOAT_LANGUAGE "HIP")
+
+    # Define `COMPAS_USE_HIP` macro so that headers can detect HIP usage
+    target_compile_definitions(${PROJECT_NAME} PUBLIC COMPAS_USE_HIP=1)
+
+    # Needed for ROCM < 6.3
+    target_compile_definitions(${PROJECT_NAME} PUBLIC HIP_ENABLE_WARP_SYNC_BUILTINS=1)
+endif()
 
 # Installation of library files
 set(COMPAS_LIB_PATH "${CMAKE_INSTALL_PREFIX}/lib")
@@ -59,14 +125,20 @@ configure_file(
     ${PROJECT_SOURCE_DIR}/julia-bindings/src/constants.h.input
     ${PROJECT_SOURCE_DIR}/julia-bindings/src/constants.h)
 
-set(CMAKE_POLICY_DEFAULT_CMP0077 NEW)
-add_subdirectory(thirdparty/Catch2)
-set(KMM_USE_CUDA ON)
 set(KMM_STATIC ON)
 add_subdirectory(thirdparty/kmm)
-target_link_libraries(${PROJECT_NAME} PUBLIC kmm)
+add_subdirectory(thirdparty/kernel_float)
 
+target_link_libraries(${PROJECT_NAME} PUBLIC kmm)
+target_link_libraries(${PROJECT_NAME} PRIVATE kernel_float)
 
 add_subdirectory(julia-bindings)
-add_subdirectory(tests)
-add_subdirectory(benchmarks)
+
+if(COMPAS_BUILD_TEST)
+    add_subdirectory(thirdparty/Catch2)
+    add_subdirectory(tests)
+endif()
+
+if(COMPAS_BUILD_BENCHMARK)
+    add_subdirectory(benchmarks)
+endif()

CompasToolkit.jl/examples/common.jl

@@ -1,9 +1,11 @@
 using BlochSimulators
 using StaticArrays
+using Statistics
 using ComputationalResources
 using ImagePhantoms
 using LinearAlgebra
 using Random
+using Printf
 
 Random.seed!(1337)
 
@@ -28,7 +30,7 @@ function generate_parameters(N)
     ρ = ComplexF32.(ImagePhantoms.shepp_logan(N, ImagePhantoms.SheppLoganEmis())') |> vec;
     T₁ = fill(0.85f0, N, N) |> vec;
     T₂ = fill(0.05f0, N, N) |> vec;
-    B₀ = repeat(1:N,1,N) .|> Float32 |> vec;
+    B₀ = zeros(Float32, N, N) |> vec; #repeat(1:N,1,N) .|> Float32 |> vec;
     B₁ = ones(Float32, N, N) |> vec;
 
     FOVˣ, FOVʸ = 25.6, 25.6;
@@ -58,6 +60,7 @@ end
 
 function generate_cartesian_trajectory(N)
     FOVˣ, FOVʸ = 25.6, 25.6;
+    FOVˣ, FOVʸ = 30.0, 28.0;
 
     nr = N # nr of readouts
     ns = N # nr of samples per readout
@@ -77,7 +80,8 @@ function generate_coils(N, ncoils)
     coil₃ = coil₁;
     coil₄ = coil₂;
 
-    return hcat(coil₁ |> vec, coil₂ |> vec, coil₃ |> vec, coil₄ |> vec) .|> Float32 
+    coils = hcat(coil₁ |> vec, coil₂ |> vec, coil₃ |> vec, coil₄ |> vec) .|> ComplexF32
+    return coils[:, 1:ncoils]
 end
 
 function generate_echos(N, sequence)
@@ -113,20 +117,45 @@ function generate_delta_echos(N, sequence)
     )
 end
 
-function print_equals_check(expected, answer)
-    atol = 1e-9
+function print_equals_check(expected, answer; atol = 1e-9)
     answer = collect(answer)
 
+    println("=== Statistics ===")
+    # Shape and length are the same, so print once
+    @printf("  Shape: %s\n", string(size(expected)))
+    @printf("  Min:   %s vs %s\n", expected[argmin(real.(expected))], answer[argmin(real.(answer))])
+    @printf("  Max:   %s vs %s\n", expected[argmax(real.(expected))], answer[argmax(real.(answer))])
+    @printf("  Mean:  %s vs %s\n", mean(expected), mean(answer))
+    @printf("  Std.:  %s vs %s\n", std(expected), std(answer))
+    @printf("  Zeros: %s vs %s\n", mean(expected .== 0), mean(answer .== 0))
+    println()
+
+    println("=== Fraction of equal values ===")
     for rtol in [0.001, 0.005, 0.01, 0.05, 0.1]
         is_equal = isapprox.(answer, expected, atol=atol, rtol=rtol)
-        println("fraction equal (atol=$(atol), rtol=$(rtol)): ", sum(is_equal) / length(answer))
+        println("  atol=$(atol), rtol=$(rtol): ", sum(is_equal) / length(answer))
     end
 
-    err = abs.(answer - expected)
-    index = argmax(err)
-    println("maximum abs error ($(index)): ", err[index], "($(answer[index]) vs $(expected[index]))")
+    err = abs.(answer .- expected)
+    idx_abs = argmax(err)
+    max_abs_err = err[idx_abs]
+
+    println("\n=== Maximum absolute error ===")
+    println("  Index:          $(Tuple(idx_abs))")
+    println("  Answer:         $(answer[idx_abs])")
+    println("  Expected:       $(expected[idx_abs])")
+    println("  Absolute Error: $(max_abs_err)")
+    println()
 
     rel_err = err ./ max.(abs.(expected), atol)
-    index = argmax(rel_err)
-    println("maximum rel error ($(index)): ", rel_err[index], "($(answer[index]) vs $(expected[index]))")
+    idx_rel = argmax(rel_err)
+    max_rel_err = rel_err[idx_rel]
+
+    println("\n=== Maximum relative error ===")
+    println("  Index:          $(Tuple(idx_rel))")
+    println("  Answer:         $(answer[idx_rel])")
+    println("  Expected:       $(expected[idx_rel])")
+    println("  Relative Error: $max_rel_err")
+    println()
+    println()
 end

CompasToolkit.jl/examples/echos.jl

@@ -1,5 +1,7 @@
 using CompasToolkit
 using ImagePhantoms
+using StructArrays
+
 include("common.jl")
 
 function make_phantom(N, coordinates)
@@ -48,6 +50,7 @@ x = LinRange(-FOVˣ/2, FOVˣ/2, N) # cm
 y = LinRange(-FOVʸ/2, FOVʸ/2, N) # cm
 
 coordinates = tuple.(x,y');
+coordinates_ref = StructArray(Coordinates(x, y, zero(x)) for (x, y) in coordinates)
 
 # Make trajectory
 # dwell time between samples within readout
@@ -77,7 +80,7 @@ parameters = vec(make_phantom(N, coordinates));
 
 # Make coil sensitivities
 ncoils = 2
-coil_sensitivities = rand((0.75:0.01:1.25), ncoils,N^2) .|> Float32
+coil_sensitivities = rand((0.75:0.01:1.25), ncoils,N^2) .|> ComplexF32
 coil_sensitivities .= 1
 coil_sensitivities = map(SVector{ncoils}, eachcol(coil_sensitivities))
 
@@ -109,14 +112,14 @@ compas_trajectory = CompasToolkit.CartesianTrajectory(
     trajectory.k_start_readout,
     trajectory.Δk_adc[1]);
 
-compas_coils = CompasToolkit.make_array(compas_context, Float32.(hcat(vec(coil₁), vec(coil₂))))
 echos = CompasToolkit.simulate_magnetization(compas_parameters, compas_sequence)
 
 # Set precision and send to gpu
 parameters          = gpu(f32(vec(parameters)))
 sequence            = gpu(f32(sequence))
 coil_sensitivities  = gpu(f32(coil_sensitivities))
 trajectory          = gpu(f32(trajectory))
+coordinates_ref     = gpu(coordinates_ref)
 
 # Simulate data
 echos_ref = simulate_magnetization(CUDALibs(), sequence, parameters)
@@ -126,6 +129,6 @@ print_equals_check(transpose(collect(echos_ref)), collect(echos))
 
 # Phase encoding
 echos = CompasToolkit.phase_encoding(echos, compas_parameters, compas_trajectory)
-phase_encoding!(echos_ref, trajectory, parameters)
+phase_encoding!(echos_ref, trajectory, coordinates_ref)
 
 print_equals_check(transpose(collect(echos_ref)), collect(echos))