diff --git a/CMakeLists.txt b/CMakeLists.txt index d486d43981..4bbab344de 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -391,9 +391,10 @@ HANDLE_SOURCES(syscheck OFF) # * FFTW (optional) Should be linked with an FFTW-like library (fftw/cufftw), # depending on whether OpenACC is enabled and which compiler is # being used. +# * LAPACK (optional) Should be linked with LAPACK. function(MFC_SETUP_TARGET) - cmake_parse_arguments(ARGS "OpenACC;MPI;SILO;HDF5;FFTW" "TARGET" "SOURCES" ${ARGN}) + cmake_parse_arguments(ARGS "OpenACC;MPI;SILO;HDF5;FFTW;LAPACK" "TARGET" "SOURCES" ${ARGN}) add_executable(${ARGS_TARGET} ${ARGS_SOURCES}) set(IPO_TARGETS ${ARGS_TARGET}) @@ -461,6 +462,11 @@ function(MFC_SETUP_TARGET) endif() endif() + if (ARGS_LAPACK) + find_package(LAPACK REQUIRED) + target_link_libraries(${a_target} PRIVATE LAPACK::LAPACK) + endif() + if (MFC_OpenACC AND ARGS_OpenACC) find_package(OpenACC) @@ -552,7 +558,7 @@ endif() if (MFC_POST_PROCESS) MFC_SETUP_TARGET(TARGET post_process SOURCES "${post_process_SRCs}" - MPI SILO HDF5 FFTW) + MPI SILO HDF5 FFTW LAPACK) # -O0 is in response to https://github.com/MFlowCode/MFC-develop/issues/95 target_compile_options(post_process PRIVATE -O0) diff --git a/src/post_process/m_lapack_example.f90 b/src/post_process/m_lapack_example.f90 new file mode 100644 index 0000000000..83aa37867b --- /dev/null +++ b/src/post_process/m_lapack_example.f90 @@ -0,0 +1,132 @@ +!> +!! @file m_lapack_example.f90 +!! @brief Contains module m_lapack_example + +!> @brief This module demonstrates the use of LAPACK in MFC post_process. +!! It provides example routines that show how to use LAPACK for +!! common linear algebra operations like solving linear systems. +module m_lapack_example + + use m_global_parameters !< Global parameters for the code + use m_mpi_proxy !< Message passing interface (MPI) module proxy + + implicit none + + private; public :: s_lapack_example_solve_linear_system, & + s_lapack_example_eigenvalues + +contains + + !> @brief Example subroutine demonstrating LAPACK usage for solving + !! a linear system Ax = b using DGESV/SGESV + !! This routine shows how to use LAPACK with MFC's precision system + impure subroutine s_lapack_example_solve_linear_system() + + ! Local variables for the linear system + integer, parameter :: n = 3 ! Size of the system + real(wp), dimension(n, n) :: A ! Coefficient matrix + real(wp), dimension(n) :: b ! Right-hand side vector + real(wp), dimension(n) :: x ! Solution vector + + ! LAPACK variables + integer, dimension(n) :: ipiv ! Pivot indices + integer :: info ! Return status + integer, parameter :: nrhs = 1 ! Number of right-hand sides + + ! Only run on the root process to avoid duplicate output + if (proc_rank /= 0) return + + ! Set up a simple 3x3 linear system: Ax = b + ! Example: + ! 2x + y + z = 8 + ! x + 3y + z = 10 + ! x + y + 4z = 16 + A(1, :) = [2.0_wp, 1.0_wp, 1.0_wp] + A(2, :) = [1.0_wp, 3.0_wp, 1.0_wp] + A(3, :) = [1.0_wp, 1.0_wp, 4.0_wp] + + b = [8.0_wp, 10.0_wp, 16.0_wp] + + print *, "LAPACK Linear System Solver Example" + print *, "Solving the system Ax = b where:" + print *, "A = [2 1 1; 1 3 1; 1 1 4]" + print *, "b = [8; 10; 16]" + + ! Copy b to x (LAPACK will overwrite the right-hand side with solution) + x = b + + ! Call appropriate LAPACK routine based on precision +#ifdef MFC_SINGLE_PRECISION + call sgesv(n, nrhs, A, n, ipiv, x, n, info) + print *, "Using LAPACK (SGESV)" +#else + call dgesv(n, nrhs, A, n, ipiv, x, n, info) + print *, "Using LAPACK (DGESV)" +#endif + + ! Check for success + if (info == 0) then + print *, "Linear system solved successfully!" + print *, "Solution: x = [", x(1), ", ", x(2), ", ", x(3), "]" + print *, "Expected: x = [1, 2, 3]" + else if (info < 0) then + print *, "LAPACK error: argument ", -info, " had an illegal value" + else + print *, "LAPACK error: matrix is singular, solution could not be computed" + end if + + print *, "End LAPACK Example" + + end subroutine s_lapack_example_solve_linear_system + + !> @brief Example subroutine demonstrating LAPACK usage for computing + !! eigenvalues of a symmetric matrix using DSYEV/SSYEV + impure subroutine s_lapack_example_eigenvalues() + + ! Local variables for eigenvalue computation + integer, parameter :: n = 3 ! Size of the matrix + real(wp), dimension(n, n) :: A ! Symmetric matrix + real(wp), dimension(n) :: w ! Eigenvalues + real(wp), dimension(3*n) :: work ! Work array + integer, parameter :: lwork = 3*n ! Size of work array + integer :: info ! Return status + character, parameter :: jobz = 'N' ! Compute eigenvalues only + character, parameter :: uplo = 'U' ! Upper triangular part of A + + ! Only run on the root process to avoid duplicate output + if (proc_rank /= 0) return + + ! Set up a simple symmetric 3x3 matrix + A(1, :) = [4.0_wp, 1.0_wp, 1.0_wp] + A(2, :) = [1.0_wp, 4.0_wp, 1.0_wp] + A(3, :) = [1.0_wp, 1.0_wp, 4.0_wp] + + print *, "LAPACK Eigenvalue Example" + print *, "Computing eigenvalues of symmetric matrix:" + print *, "A = [4 1 1; 1 4 1; 1 1 4]" + + ! Call appropriate LAPACK routine based on precision +#ifdef MFC_SINGLE_PRECISION + call ssyev(jobz, uplo, n, A, n, w, work, lwork, info) + print *, "Using LAPACK (SSYEV)" +#else + call dsyev(jobz, uplo, n, A, n, w, work, lwork, info) + print *, "Using LAPACK (DSYEV)" +#endif + + ! Check for success + if (info == 0) then + print *, "Eigenvalues computed successfully!" + print *, "Eigenvalues: [", w(1), ", ", w(2), ", ", w(3), "]" + print *, "Expected: [2, 5, 5] (approximately)" + else if (info < 0) then + print *, "LAPACK error: argument ", -info, " had an illegal value" + else + print *, "LAPACK error: algorithm failed to converge" + end if + + print *, "End LAPACK Eigenvalue Example" + + end subroutine s_lapack_example_eigenvalues + +end module m_lapack_example diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index b733f43dd9..12162d1516 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -45,6 +45,8 @@ module m_start_up use m_chemistry + use m_lapack_example + implicit none contains @@ -714,6 +716,11 @@ impure subroutine s_initialize_mpi_domain ! leads to the termination of the post-process. if (proc_rank == 0) then call s_assign_default_values_to_user_inputs() + + ! Run LAPACK examples before reading input file + call s_lapack_example_solve_linear_system() + call s_lapack_example_eigenvalues() + call s_read_input_file() call s_check_input_file() diff --git a/toolchain/cmake/regular/FindLAPACK.cmake b/toolchain/cmake/regular/FindLAPACK.cmake new file mode 100644 index 0000000000..6691eb40bb --- /dev/null +++ b/toolchain/cmake/regular/FindLAPACK.cmake @@ -0,0 +1,60 @@ +# Attempt to find LAPACK (Linear Algebra PACKage) +# URL: https://www.netlib.org/lapack/ +# DOCS: https://cmake.org/cmake/help/latest/command/find_library.html +# https://cmake.org/cmake/help/latest/module/FindPackageHandleStandardArgs.html + +include(FindPackageHandleStandardArgs) + +# Special handling for Cray systems which have optimized math libraries +if (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray") + # On Cray systems, LAPACK is typically provided by the cray-libsci package + find_library(LAPACK_LIBRARY + NAMES sci_cray sci_gnu sci_intel sci_pgi sci + NAMES_PER_DIR + ) + set(BLAS_LIBRARY "") # BLAS is included in the sci library +else() + # Find LAPACK library for other compilers + find_library(LAPACK_LIBRARY + NAMES lapack + PATH_SUFFIXES lapack + NAMES_PER_DIR + ) + + # Find BLAS library (required by LAPACK) + find_library(BLAS_LIBRARY + NAMES blas openblas + PATH_SUFFIXES blas + NAMES_PER_DIR + ) + + # Some LAPACK implementations include BLAS + if (NOT BLAS_LIBRARY) + set(BLAS_LIBRARY "") + endif() +endif() + +FIND_PACKAGE_HANDLE_STANDARD_ARGS( + LAPACK + REQUIRED_VARS + LAPACK_LIBRARY +) + +if (LAPACK_FOUND AND NOT TARGET LAPACK::LAPACK) + set(LAPACK_LIBRARIES "${LAPACK_LIBRARY}") + if (BLAS_LIBRARY) + list(APPEND LAPACK_LIBRARIES "${BLAS_LIBRARY}") + endif() + + add_library(LAPACK::LAPACK INTERFACE IMPORTED) + + set_target_properties(LAPACK::LAPACK PROPERTIES + INTERFACE_LINK_LIBRARIES "${LAPACK_LIBRARIES}" + ) + + # Add math library for linking (commonly needed) + if (UNIX AND NOT APPLE) + set_property(TARGET LAPACK::LAPACK APPEND PROPERTY + INTERFACE_LINK_LIBRARIES "m") + endif() +endif() \ No newline at end of file diff --git a/toolchain/dependencies/CMakeLists.txt b/toolchain/dependencies/CMakeLists.txt index 116e288711..428bcb7204 100644 --- a/toolchain/dependencies/CMakeLists.txt +++ b/toolchain/dependencies/CMakeLists.txt @@ -21,6 +21,7 @@ include(ExternalProject) option(MFC_FFTW "Build the FFTW3 dependency" OFF) option(MFC_HDF5 "Build the HDF5 dependency" OFF) option(MFC_SILO "Build the SILO dependency" OFF) +option(MFC_LAPACK "Build the LAPACK dependency" OFF) option(MFC_HIPFORT "Build the HIPFORT dependency" OFF) @@ -98,6 +99,32 @@ if (MFC_SILO) endif() endif() +# LAPACK +if (MFC_LAPACK) + find_package(LAPACK QUIET) + if (LAPACK_FOUND) + message(STATUS "LAPACK found.") + add_custom_target(lapack) + else() + # Use reference LAPACK from netlib + find_package(Git REQUIRED) + + ExternalProject_Add(lapack + GIT_REPOSITORY "https://github.com/Reference-LAPACK/lapack" + GIT_TAG v3.12.0 + GIT_SHALLOW ON + GIT_PROGRESS ON + CMAKE_ARGS -DBUILD_SHARED_LIBS=OFF + -DBUILD_TESTING=OFF + -DCBLAS=OFF + -DLAPACKE=OFF + -DBUILD_DEPRECATED=OFF + "-DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX}" + "-DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}" + ) + endif() +endif() + # HIPFORT if (MFC_HIPFORT) if (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray") diff --git a/toolchain/mfc/build.py b/toolchain/mfc/build.py index 2de738986d..b28be82821 100644 --- a/toolchain/mfc/build.py +++ b/toolchain/mfc/build.py @@ -184,14 +184,15 @@ def install(self, case: input.MFCInputFile): FFTW = MFCTarget('fftw', ['-DMFC_FFTW=ON'], True, False, False, MFCTarget.Dependencies([], [], []), -1) HDF5 = MFCTarget('hdf5', ['-DMFC_HDF5=ON'], True, False, False, MFCTarget.Dependencies([], [], []), -1) SILO = MFCTarget('silo', ['-DMFC_SILO=ON'], True, False, False, MFCTarget.Dependencies([HDF5], [], []), -1) +LAPACK = MFCTarget('lapack', ['-DMFC_LAPACK=ON'], True, False, False, MFCTarget.Dependencies([], [], []), -1) HIPFORT = MFCTarget('hipfort', ['-DMFC_HIPFORT=ON'], True, False, False, MFCTarget.Dependencies([], [], []), -1) PRE_PROCESS = MFCTarget('pre_process', ['-DMFC_PRE_PROCESS=ON'], False, True, False, MFCTarget.Dependencies([], [], []), 0) SIMULATION = MFCTarget('simulation', ['-DMFC_SIMULATION=ON'], False, True, False, MFCTarget.Dependencies([], [FFTW], [HIPFORT]), 1) -POST_PROCESS = MFCTarget('post_process', ['-DMFC_POST_PROCESS=ON'], False, True, False, MFCTarget.Dependencies([FFTW, HDF5, SILO], [], []), 2) +POST_PROCESS = MFCTarget('post_process', ['-DMFC_POST_PROCESS=ON'], False, True, False, MFCTarget.Dependencies([FFTW, HDF5, SILO, LAPACK], [], []), 2) SYSCHECK = MFCTarget('syscheck', ['-DMFC_SYSCHECK=ON'], False, False, True, MFCTarget.Dependencies([], [], [HIPFORT]), -1) DOCUMENTATION = MFCTarget('documentation', ['-DMFC_DOCUMENTATION=ON'], False, False, False, MFCTarget.Dependencies([], [], []), -1) -TARGETS = { FFTW, HDF5, SILO, HIPFORT, PRE_PROCESS, SIMULATION, POST_PROCESS, SYSCHECK, DOCUMENTATION } +TARGETS = { FFTW, HDF5, SILO, LAPACK, HIPFORT, PRE_PROCESS, SIMULATION, POST_PROCESS, SYSCHECK, DOCUMENTATION } DEFAULT_TARGETS = { target for target in TARGETS if target.isDefault } REQUIRED_TARGETS = { target for target in TARGETS if target.isRequired }