Replace Numerial Recipes LU-decomp routines with LAPACK ones

.github/workflows/Intel.yml

@@ -37,7 +37,7 @@ jobs:
         rm GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB
         echo "deb https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list
         sudo apt-get update
-        sudo apt-get install intel-oneapi-openmp intel-oneapi-compiler-fortran-2023.2.1 intel-oneapi-compiler-dpcpp-cpp-and-cpp-classic-2023.2.1
+        sudo apt-get install intel-oneapi-openmp intel-oneapi-compiler-fortran-2023.2.1 intel-oneapi-compiler-dpcpp-cpp-and-cpp-classic-2023.2.1 intel-oneapi-mkl-devel-2023.2.0
         echo "source /opt/intel/oneapi/setvars.sh" >> ~/.bash_profile
 
     - name: checkout

.github/workflows/Linux.yml

@@ -33,6 +33,7 @@ jobs:
 
     - name: build
       run: |
+        sudo apt install libopenblas-dev
         cd ip
         mkdir build 
         cd build

.github/workflows/Spack.yml

@@ -33,6 +33,7 @@ jobs:
 
     - name: spack-build-and-test
       run: |
+        sudo apt install libopenblas-dev
         git clone -c feature.manyFiles=true https://github.com/spack/spack
         . spack/share/spack/setup-env.sh
         spack env create ip-env
@@ -42,7 +43,7 @@ jobs:
         spack add ip@develop%gcc@11 ${{ matrix.variants }} target=x86_64
         precision=$(echo ${{ matrix.variants }} | grep -oP " precision=\K[4d8]")
         if [ "$precision" == "d" ]; then spack add grib-util@develop ; fi
-        spack external find cmake gmake
+        spack external find cmake gmake openblas
         spack concretize
         # Run installation and run CTest suite
         spack install --verbose --fail-fast --test root

.github/workflows/developer.yml

@@ -25,7 +25,7 @@ jobs:
     - name: Install Dependencies
       run: |
         sudo apt-get update
-        sudo apt-get install doxygen
+        sudo apt-get install doxygen libopenblas-dev
         python3 -m pip install gcovr
 
     - name: checkout

CMakeLists.txt

@@ -49,6 +49,11 @@ if(OPENMP)
   find_package(OpenMP REQUIRED COMPONENTS Fortran)
 endif()
 
+if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$")
+  set(BLA_VENDOR Intel10_64lp_seq)
+endif()
+find_package(LAPACK)
+
 # Set compiler flags.
 if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$")
   set(CMAKE_Fortran_FLAGS "-g -traceback -assume byterecl -fp-model strict -fpp -auto ${CMAKE_Fortran_FLAGS}")

spack/package.py

@@ -65,6 +65,7 @@ class Ip(CMakePackage):
     depends_on("sp precision=4", when="@4.1:4 precision=4")
     depends_on("sp precision=d", when="@4.1:4 precision=d")
     depends_on("sp precision=8", when="@4.1:4 precision=8")
+    depends_on("lapack", when="@develop")
 
     def cmake_args(self):
         args = [

src/CMakeLists.txt

@@ -14,7 +14,7 @@ ip_mercator_grid_mod.F90 ip_polar_stereo_grid_mod.F90
 ip_rot_equid_cylind_egrid_mod.F90 ip_rot_equid_cylind_grid_mod.F90
 ip_constants_mod.F90 ip_grids_mod.F90 ip_grid_factory_mod.F90
 ip_interpolators_mod.F90 earth_radius_mod.F90 polfix_mod.F90
-fftpack.F lapack_gen.F ncpus.F spanaly.f spdz2uv.f speps.f spfft1.f spffte.f
+fftpack.F ncpus.F spanaly.f spdz2uv.f speps.f spfft1.f spffte.f
 spfftpt.f splaplac.f splat.F splegend.f sppad.f spsynth.f sptezd.f sptez.f
 sptezmd.f sptezm.f sptezmv.f sptezv.f sptgpm.f sptgpmv.f sptgps.f sptgpsv.f
 sptgpt.f sptgptv.f sptrand.f sptran.f sptranf0.f sptranf1.f sptranf.f sptranfv.f
@@ -77,6 +77,7 @@ foreach(kind ${kinds})
   if(OpenMP_Fortran_FOUND)
     target_link_libraries(${lib_name} PUBLIC OpenMP::OpenMP_Fortran)
   endif()
+  target_link_libraries(${lib_name} PUBLIC LAPACK::LAPACK)
 
   list(APPEND LIB_TARGETS ${lib_name})
 

src/splat.F

@@ -64,7 +64,7 @@ SUBROUTINE SPLAT(IDRT,JMAX,SLAT,WLAT)
      $ 146.870307625, 150.011882457, 153.153458019, 156.295034268 /
       REAL:: DLT,D1=1.
       REAL AWORK((JMAX+1)/2,((JMAX+1)/2)),BWORK(((JMAX+1)/2))
-      INTEGER:: JHE,JHO,J0=0
+      INTEGER:: JHE,JHO,J0=0, INFO
       INTEGER IPVT((JMAX+1)/2)
       PARAMETER(PI=3.14159265358979,C=(1.-(2./PI)**2)*0.25)
 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@@ -133,8 +133,14 @@ SUBROUTINE SPLAT(IDRT,JMAX,SLAT,WLAT)
           BWORK(JS)=-D1/(4*(JS-1)**2-1)
         ENDDO
 
-        call ludcmp(awork,jho,jhe,ipvt)
-        call lubksb(awork,jho,jhe,ipvt,bwork)
+        ! Call LAPACK routines
+#if (LSIZE==4)
+        CALL SGETRF(JHO, JHO, AWORK, JHE, IPVT, INFO)
+        CALL SGETRS('N', JHO, 1, AWORK, JHE, IPVT, BWORK, JHO, INFO)
+#else
+        CALL DGETRF(JHO, JHO, AWORK, JHE, IPVT, INFO)
+        CALL DGETRS('N', JHO, 1, AWORK, JHE, IPVT, BWORK, JHO, INFO)
+#endif
 
         WLAT(1)=0.
         DO J=1,JHO
@@ -170,8 +176,14 @@ SUBROUTINE SPLAT(IDRT,JMAX,SLAT,WLAT)
           BWORK(JS)=-D1/(4*(JS-1)**2-1)
         ENDDO
 
-        call ludcmp(awork,jho,jhe,ipvt,d)
-        call lubksb(awork,jho,jhe,ipvt,bwork)
+        ! Call LAPACK routines
+#if (LSIZE==4)
+        CALL SGETRF(JHO, JHO, AWORK, JHE, IPVT, INFO)
+        CALL SGETRS('N', JHO, 1, AWORK, JHE, IPVT, BWORK, JHO, INFO)
+#else
+        CALL DGETRF(JHO, JHO, AWORK, JHE, IPVT, INFO)
+        CALL DGETRS('N', JHO, 1, AWORK, JHE, IPVT, BWORK, JHO, INFO)
+#endif
 
         WLAT(1)=0.
         DO J=1,JHO