Feature: Fortran interface for lapack dense solver (#353)

Added Fortran support for the LAPACK dense SUNLinearSolver implementation.

---------

Co-authored-by: Daniel R. Reynolds <reynolds@smu.edu>

CHANGELOG.md

@@ -27,6 +27,8 @@ Fixed scaling bug in `SUNMatScaleAddI_Sparse` for non-square matrices.
 Fixed missing soversions in some `SUNLinearSolver` and `SUNNonlinearSolver`
 CMake targets.
 
+Added Fortran support for the LAPACK dense `SUNLinearSolver` implementations.
+
 Added the fourth order ERK method `ARKODE_SOFRONIOU_SPALETTA_5_3_4`.
 
 ## Changes to SUNDIALS in release 6.6.2

doc/arkode/guide/source/Introduction.rst

@@ -155,6 +155,8 @@ CMake targets.
 
 Added the fourth order ERK method ``ARKODE_SOFRONIOU_SPALETTA_5_3_4``.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v5.6.2
 -----------------
 

doc/cvode/guide/source/Introduction.rst

@@ -126,6 +126,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v6.6.2
 -----------------
 

doc/cvodes/guide/source/Introduction.rst

@@ -129,6 +129,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v6.6.2
 -----------------
 

doc/ida/guide/source/Introduction.rst

@@ -87,6 +87,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v6.6.2
 -----------------
 

doc/idas/guide/source/Introduction.rst

@@ -104,6 +104,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v5.6.2
 -----------------
 

doc/kinsol/guide/source/Introduction.rst

@@ -99,6 +99,8 @@ Fixed scaling bug in ``SUNMatScaleAddI_Sparse`` for non-square matrices.
 Fixed missing soversions in some ``SUNLinearSolver`` and ``SUNNonlinearSolver``
 CMake targets.
 
+Added Fortran support for the LAPACK  dense ``SUNLinearSolver`` implementation.
+
 Changes in v6.6.2
 -----------------
 

examples/sunlinsol/lapackdense/CMakeLists.txt

@@ -30,6 +30,16 @@ set(sunlinsol_lapackdense_dependencies
   test_sunlinsol
 )
 
+# If building F2003 tests
+if (BUILD_FORTRAN_MODULE_INTERFACE)
+  set(sunlinsol_lapackdense_fortran_examples
+    "test_fsunlinsol_lapackdense_mod\;10 0 0\;"
+    "test_fsunlinsol_lapackdense_mod\;100 0 0\;"
+    "test_fsunlinsol_lapackdense_mod\;500 0 0\;"
+    "test_fsunlinsol_lapackdense_mod\;1000 0 0\;"
+    )
+endif()
+
 # Add source directory to include directories
 include_directories(. ..)
 
@@ -82,6 +92,52 @@ foreach(example_tuple ${sunlinsol_lapackdense_examples})
 
 endforeach(example_tuple ${sunlinsol_lapackdense_examples})
 
+# Add the build and install targets for each example
+foreach(example_tuple ${sunlinsol_lapackdense_fortran_examples})
+
+  # parse the example tuple
+  list(GET example_tuple 0 example)
+  list(GET example_tuple 1 example_args)
+  list(GET example_tuple 2 example_type)
+
+  # check if this example has already been added, only need to add
+  # example source files once for testing with different inputs
+  if(NOT TARGET ${example})
+    # example source files
+    add_executable(${example} ${example}.f90
+      ${SUNDIALS_SOURCE_DIR}/examples/utilities/test_utilities.f90
+      ${SUNDIALS_SOURCE_DIR}/examples/sunlinsol/test_sunlinsol.f90)
+
+    # folder to organize targets in an IDE
+    set_target_properties(${example} PROPERTIES FOLDER "Examples")
+
+    # set fortran module directory to avoid name collisions
+    set_target_properties(${example} PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+
+    # libraries to link against
+    target_link_libraries(${example}
+      sundials_nvecserial
+      sundials_fnvecserial_mod
+      sundials_sunlinsollapackdense
+      sundials_fsunlinsollapackdense_mod
+      ${EXE_EXTRA_LINK_LIBS})
+
+    # check if example args are provided and set the test name
+    if("${example_args}" STREQUAL "")
+      set(test_name ${example})
+    else()
+      string(REGEX REPLACE " " "_" test_name ${example}_${example_args})
+    endif()
+
+    # add example to regression tests
+    sundials_add_test(${test_name} ${example}
+      TEST_ARGS ${example_args}
+      EXAMPLE_TYPE ${example_type}
+      NODIFF)
+  endif()
+
+endforeach(example_tuple ${sunlinsol_lapackdense_fortran_examples})
+
 
 if(EXAMPLES_INSTALL)
 

examples/sunlinsol/lapackdense/test_fsunlinsol_lapackdense_mod.f90

@@ -0,0 +1,185 @@
+! -----------------------------------------------------------------
+! Programmer(s): Daniel R. Reynolds @ SMU
+! -----------------------------------------------------------------
+! SUNDIALS Copyright Start
+! Copyright (c) 2002-2023, Lawrence Livermore National Security
+! and Southern Methodist University.
+! All rights reserved.
+!
+! See the top-level LICENSE and NOTICE files for details.
+!
+! SPDX-License-Identifier: BSD-3-Clause
+! SUNDIALS Copyright End
+! -----------------------------------------------------------------
+! This file tests the Fortran 2003 interface to the SUNDIALS
+! LAPACK-Dense SUNLinearSolver implementation.
+! -----------------------------------------------------------------
+
+module test_fsunlinsol_lapackdense
+  use, intrinsic :: iso_c_binding
+  use test_utilities
+  implicit none
+
+  integer(C_LONG), private, parameter :: N = 100
+
+contains
+
+  integer(C_INT) function unit_tests() result(fails)
+    use, intrinsic :: iso_c_binding
+    use fsundials_nvector_mod
+    use fsundials_matrix_mod
+    use fsundials_linearsolver_mod
+    use fnvector_serial_mod
+    use fsunmatrix_dense_mod
+    use fsunlinsol_lapackdense_mod
+    use test_sunlinsol
+
+    implicit none
+
+    type(SUNLinearSolver), pointer :: LS                ! test linear solver
+    type(SUNMatrix),       pointer :: A, I              ! test matrices
+    type(N_Vector),        pointer :: x, b              ! test vectors
+    real(C_DOUBLE),        pointer :: colj(:), colIj(:) ! matrix column data
+    real(C_DOUBLE),        pointer :: xdata(:)          ! x vector data
+    real(C_DOUBLE)                 :: tmpr              ! temporary real value
+    integer(C_LONG)                :: j, k
+    integer(C_INT)                 :: tmp
+
+    fails = 0
+
+    A => FSUNDenseMatrix(N, N, sunctx)
+    I => FSUNDenseMatrix(N, N, sunctx)
+    x => FN_VNew_Serial(N, sunctx)
+    b => FN_VNew_Serial(N, sunctx)
+
+    ! fill A matrix with uniform random data in [0, 1/N)
+    do j=1, N
+      colj => FSUNDenseMatrix_Column(A, j-1)
+      do k=1, N
+        call random_number(tmpr)
+        colj(k) = tmpr / N
+      end do
+    end do
+
+    ! create anti-identity matrix
+    j = N
+    do k=1, N
+      colj => FSUNDenseMatrix_Column(I, j-1)
+      colj(k) = ONE
+      j = j-1
+    end do
+
+    ! add anti-identity to ensure the solver needs to do row-swapping
+    do k=1, N
+      do j=1, N
+        colj => FSUNDenseMatrix_Column(A, j-1)
+        colIj => FSUNDenseMatrix_Column(I, j-1)
+        colj(k) = colj(k) + colIj(k)
+      end do
+    end do
+
+    ! fill x vector with uniform random data in [0, 1)
+    xdata => FN_VGetArrayPointer(x)
+    do j=1, N
+      call random_number(tmpr)
+      xdata(j) = tmpr
+    end do
+
+    ! create RHS vector for linear solve
+    fails = FSUNMatMatvec(A, x, b)
+    if (fails /= 0) then
+      call FSUNMatDestroy(A)
+      call FSUNMatDestroy(I)
+      call FN_VDestroy(x)
+      call FN_VDestroy(b)
+      return
+    end if
+
+    ! create dense linear solver
+    LS => FSUNLinSol_LapackDense(x, A, sunctx)
+
+    ! run tests
+    fails = fails + Test_FSUNLinSolInitialize(LS, 0)
+    fails = fails + Test_FSUNLinSolSetup(LS, A, 0)
+    fails = fails + Test_FSUNLinSolSolve(LS, A, x, b, 100*UNIT_ROUNDOFF, 0)
+
+    fails = fails + Test_FSUNLinSolGetType(LS, SUNLINEARSOLVER_DIRECT, 0)
+    fails = fails + Test_FSUNLinSolLastFlag(LS, 0)
+    fails = fails + Test_FSUNLinSolSpace(LS, 0)
+
+    ! cleanup
+    tmp = FSUNLinSolFree(LS)
+    call FSUNMatDestroy(A)
+    call FSUNMatDestroy(I)
+    call FN_VDestroy(x)
+    call FN_VDestroy(b)
+
+  end function unit_tests
+
+end module
+
+integer(C_INT) function check_vector(X, Y, tol) result(failure)
+  use, intrinsic :: iso_c_binding
+  use fsundials_nvector_mod
+  use test_utilities
+
+  implicit none
+  type(N_Vector)  :: x, y
+  real(C_DOUBLE)  :: tol, maxerr
+  integer(C_LONG) :: i, xlen, ylen
+  real(C_DOUBLE), pointer :: xdata(:), ydata(:)
+
+  failure = 0
+
+  xdata => FN_VGetArrayPointer(x)
+  ydata => FN_VGetArrayPointer(y)
+
+  xlen = FN_VGetLength(x)
+  ylen = FN_VGetLength(y)
+
+  if (xlen /= ylen) then
+    print *, 'FAIL: check_vector: different data array lengths'
+    failure = 1
+    return
+  end if
+
+  do i = 1, xlen
+    failure = failure + FNEQTOL(xdata(i), ydata(i), tol)
+  end do
+
+  if (failure > 0) then
+    maxerr = ZERO
+    do i = 1, xlen
+      maxerr = max(abs(xdata(i)-ydata(i))/abs(ydata(i)), maxerr)
+    end do
+    write(*,'(A,E14.7,A,E14.7,A)') &
+      "FAIL: check_vector failure: maxerr = ", maxerr, "  (tol = ", tol, ")"
+  end if
+
+end function check_vector
+
+program main
+  !======== Inclusions ==========
+  use, intrinsic :: iso_c_binding
+  use test_fsunlinsol_lapackdense
+
+  !======== Declarations ========
+  implicit none
+  integer(C_INT) :: fails = 0
+
+  !============== Introduction =============
+  print *, 'LAPACK-Dense SUNLinearSolver Fortran 2003 interface test'
+
+  call Test_Init(c_null_ptr)
+
+  fails = unit_tests()
+  if (fails /= 0) then
+    print *, 'FAILURE: n unit tests failed'
+    stop 1
+  else
+    print *,'SUCCESS: all unit tests passed'
+  end if
+
+  call Test_Finalize()
+
+end program main

scripts/shared

@@ -490,6 +490,9 @@ echo "   --- Add sunlinsol/lapackdense module to $tarfile"
 
 $tar $tarfile $distrobase/src/sunlinsol/lapackdense/CMakeLists.txt
 $tar $tarfile $distrobase/src/sunlinsol/lapackdense/sunlinsol_lapackdense.c
+$tar $tarfile $distrobase/src/sunlinsol/lapackdense/fmod/CMakeLists.txt
+$tar $tarfile $distrobase/src/sunlinsol/lapackdense/fmod/fsunlinsol_lapackdense_mod.c
+$tar $tarfile $distrobase/src/sunlinsol/lapackdense/fmod/fsunlinsol_lapackdense_mod.f90
 
 echo "   --- Add sunlinsol/magmadense module to $tarfile"
 
@@ -713,6 +716,7 @@ $tar $tarfile $distrobase/examples/sunlinsol/kokkos/CMakeLists.txt
 $tar $tarfile $distrobase/examples/sunlinsol/lapackband/test_sunlinsol_lapackband.c
 $tar $tarfile $distrobase/examples/sunlinsol/lapackband/CMakeLists.txt
 $tar $tarfile $distrobase/examples/sunlinsol/lapackdense/test_sunlinsol_lapackdense.c
+$tar $tarfile $distrobase/examples/sunlinsol/lapackdense/test_fsunlinsol_lapackdense_mod.f90
 $tar $tarfile $distrobase/examples/sunlinsol/lapackdense/CMakeLists.txt
 
 $tar $tarfile $distrobase/examples/sunlinsol/magmadense/test_sunlinsol_magmadense.cpp

src/sunlinsol/lapackdense/CMakeLists.txt

@@ -37,3 +37,8 @@ sundials_add_library(sundials_sunlinsollapackdense
 )
 
 message(STATUS "Added SUNLINSOL_LAPACKDENSE module")
+
+# Add module if F2003 interface is enabled
+if(BUILD_FORTRAN_MODULE_INTERFACE)
+  add_subdirectory(fmod)
+endif()

src/sunlinsol/lapackdense/fmod/CMakeLists.txt

@@ -0,0 +1,32 @@
+# ----------------------------------------------------------------------
+# Programmer(s): Cody J. Balos @ LLNL
+# ----------------------------------------------------------------------
+# SUNDIALS Copyright Start
+# Copyright (c) 2002-2023, Lawrence Livermore National Security
+# and Southern Methodist University.
+# All rights reserved.
+#
+# See the top-level LICENSE and NOTICE files for details.
+#
+# SPDX-License-Identifier: BSD-3-Clause
+# SUNDIALS Copyright End
+# ----------------------------------------------------------------------
+# CMakeLists.txt file for the F2003 LAPACK dense SUNLinearSolver 
+# object library
+# ----------------------------------------------------------------------
+
+sundials_add_f2003_library(sundials_fsunlinsollapackdense_mod
+  SOURCES
+    fsunlinsol_lapackdense_mod.f90 fsunlinsol_lapackdense_mod.c
+  OBJECT_LIBRARIES
+    sundials_fgeneric_mod_obj
+  LINK_LIBRARIES
+    PUBLIC sundials_fsunmatrixdense_mod
+  OUTPUT_NAME
+    sundials_fsunlinsollapackdense_mod
+  VERSION
+    ${sunlinsollib_VERSION}
+  SOVERSION
+    ${sunlinsollib_SOVERSION}
+)
+message(STATUS "Added SUNLINSOL_LAPACKDENSE F2003 interface")