Example of using LAPACK in post_process

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)

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

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()
 

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() 

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")

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 }