Support of BLAS helpers

.github/workflows/CI.yml

@@ -26,9 +26,10 @@ jobs:
 #          - {compiler: gcc, version: 11}
 #          - {compiler: gcc, version: 12}
 #          - {compiler: gcc, version: 13}
-          - {compiler: intel, version: '2024.2', mkl: '2024.2'}
-          - {compiler: intel, version: '2024.1', mkl: '2024.1'}
-          - {compiler: intel-classic, version: '2021.9', mkl: '2021.4.0'}
+          - {compiler: intel, version: '2024.2', mkl: '2024.2', mkl_sparse: '1'}
+          - {compiler: intel, version: '2024.1', mkl: '2024.1', mkl_sparse: '1'}
+          - {compiler: intel, version: latest, mkl: '2026.1', mkl_sparse: '0'}
+          - {compiler: intel-classic, version: '2021.9', mkl: '2021.4.0', mkl_sparse: '1'}
         build: [cmake]
     env:
       BUILD_DIR: ${{ matrix.build == 'cmake' && 'build' || '.' }}
@@ -46,7 +47,7 @@ jobs:
       run: pip install --upgrade fypp ninja
 
     - name: Setup Fortran compiler
-      uses: fortran-lang/setup-fortran@v1.6.1
+      uses: fortran-lang/setup-fortran@v1.9.2
       id: setup-fortran
       with:
         compiler: ${{ matrix.toolchain.compiler }}
@@ -65,6 +66,7 @@ jobs:
       if: ${{ contains(matrix.build, 'cmake') }}
       run: >-
         cmake -Wdev
+        -DCMAKE_MKL_SPARSE=${{ matrix.toolchain.mkl_sparse }}
         -DCMAKE_SPAINV=0
         -DCMAKE_BUILD_TYPE=Debug
         -DCMAKE_INSTALL_PREFIX=$PWD/_dist

CMakeLists.txt

@@ -89,6 +89,14 @@ endif()
 #endif()
 message(STATUS "Support of MKL PARDISO: ${WITH_MKLPARDISO}")
 
+#MKL SPARSE
+if(NOT DEFINED CMAKE_MKL_SPARSE)
+  set(WITH_MKL_SPARSE 1 CACHE STRING "Support of MKL SPARSE")
+else()
+  set(WITH_MKL_SPARSE ${CMAKE_MKL_SPARSE} CACHE STRING "Support of MKL SPARSE")
+endif()
+
+
 #SPAINV
 if(NOT DEFINED CMAKE_SPAINV)
   set(WITH_SPAINV 0 CACHE STRING "Support of SPAINV")
@@ -145,6 +153,7 @@ target_compile_definitions(
   PRIVATE
   "-D_DP=1"
   "-D_METIS=${WITH_METIS}"
+  "-D_MKL_SPARSE=${WITH_MKL_SPARSE}"
   "-D_PARDISO=${WITH_MKLPARDISO}"
   "-D_SPAINV=${WITH_SPAINV}"
   "-D_VERBOSE=${WITH_VERBOSE}"

Makefile

@@ -3,6 +3,7 @@
 DEBUGENABLE=0
 DPENABLE=1
 METISENABLE=0
+MKLSPARSEENABLE=1
 PARDISOENABLE=1
 SPAINVENABLE=1
 VERBOSE=0
@@ -63,6 +64,12 @@ else
  METIS = 0
 endif
 
+ifeq ($(MKLSPARSEENABLE), 1)
+ MKL_SPARSE=1
+else
+ MKL_SPARSE=0
+endif
+
 
 ifeq ($(DPENABLE),0)
  DP=0
@@ -83,7 +90,7 @@ else
  SPAINV=0
 endif
 
-FFLAGS += -D_DP=$(DP) -D_METIS=$(METIS) -D_PARDISO=$(PARDISO) -D_SPAINV=$(SPAINV) -D_VERBOSE=$(VERBOSE)
+FFLAGS += -D_DP=$(DP) -D_METIS=$(METIS) -D_MKL_SPARSE=$(MKL_SPARSE) -D_PARDISO=$(PARDISO) -D_SPAINV=$(SPAINV) -D_VERBOSE=$(VERBOSE)
 
 FYPPFLAGS =
 

src/CMakeLists.txt

@@ -26,6 +26,7 @@ list(
   "${dir}/modsparse_crs.f90"
   "${dir}/modsparse.f90"
   "${dir}/modsparse_gen.f90"
+  "${dir}/modsparse_helpers.f90"
   "${dir}/modsparse_inv.f90"
   "${dir}/modsparse_inv_int.f90"
   "${dir}/modsparse_ll.f90"

src/Makefile

@@ -20,6 +20,7 @@ OBJ = modcommon.o \
       modsparse_inv_int.o \
       modsparse_inv.o  smbfct.o \
       $(OBJPARDISO) $(OBJMETIS) $(OBJSPAINV) \
+      modsparse_helpers.o \
       modsparse_mkl.o modsparse_hash.o modsparse.o modsparse_gen.o modsparse_coo.o \
       modsparse_crs64.o \
       modsparse_crs.o modsparse_ll.o modsparse_metisgraph.o
@@ -36,7 +37,7 @@ libsparse.a: $(OBJ)
 
 #dependencies
 modsparse_coo.o: modsparse.o modsparse_hash.o
-modsparse_crs.o: modsparse_mkl.o modsparse.o $(OBJSPAINV) $(OBJPARDISO) $(OBJMETIS)
+modsparse_crs.o: modsparse_mkl.o modsparse.o modsparse_helpers.o $(OBJSPAINV) $(OBJPARDISO) $(OBJMETIS)
 modsparse_crs64.o: modsparse.o $(OBJPARDISO)
 ifeq ($(SPAINVENABLE),1)
 modsparse_inv.o: modspainv.o smbfct.o modsparse_inv_int.o

src/modsparse_crs.f90

@@ -1,10 +1,14 @@
 submodule (modsparse) modsparse_crs
  use modrandom, only: setseed, snorm=>rand_stdnormal
- use modsparse_mkl, only: pardisoinit, pardiso &
-                          , mkl_scsrmv, mkl_dcsrmv &
+ use modsparse_mkl, only: pardisoinit, pardiso
+#if (_MKL_SPARSE==1)
+ use modsparse_mkl, only: mkl_scsrmv, mkl_dcsrmv &
                           , mkl_scsrmm, mkl_dcsrmm &
                           , mkl_scsrtrsv, mkl_dcsrtrsv &
                           , mkl_scsrsymv, mkl_dcsrsymv
+#else
+ use modsparse_helpers, only: csrsymv, csrmm, csrmv, csrtrsv
+#endif
  use modsparse_inv, only: get_chol, get_ichol, get_spainv
 #if (_PARDISO==1)
  use modvariablepardiso, only: checkpardiso, pardiso_variable
@@ -251,10 +255,14 @@ module subroutine harville_crs(sparse, ngibbs, nburn, diaginv, seed)
  endif
 
  do i = 1, ngibbs
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
   call mkl_scsrsymv('U', n, a, sparse%ia, sparse%ja, xt, x)
 #else
   call mkl_dcsrsymv('U', n, a, sparse%ia, sparse%ja, xt, x)
+#endif
+#else
+  call csrsymv('U', n, a, sparse%ia, sparse%ja, xt, x)
 #endif
   x = -1 * x
   do j=1,n
@@ -339,6 +347,7 @@ module subroutine multgenv_csr(sparse,alpha,trans,x,val,y)
 
  matdescra(4)='F'
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrmv(trans,sparse%dim1,sparse%dim2,alpha,matdescra&
                   ,sparse%a,sparse%ja,sparse%ia(1:sparse%dim1),sparse%ia(2:sparse%dim1+1)&
@@ -348,6 +357,11 @@ module subroutine multgenv_csr(sparse,alpha,trans,x,val,y)
                   ,sparse%a,sparse%ja,sparse%ia(1:sparse%dim1),sparse%ia(2:sparse%dim1+1)&
                   ,x,val,y)
 #endif
+#else
+ call csrmv(trans,sparse%dim1,sparse%dim2,alpha,matdescra&
+                  ,sparse%a,sparse%ja,sparse%ia(1:sparse%dim1),sparse%ia(2:sparse%dim1+1)&
+                  ,x,val,y)
+#endif
 
 end subroutine
 
@@ -382,6 +396,7 @@ module subroutine multgenm_csr(sparse,alpha,trans,x,val,y)
 
  matdescra(4)='F'
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrmm(trans,sparse%dim1,size(y,2),sparse%dim2,&
                  alpha,matdescra,sparse%a,sparse%ja,sparse%ia(1:sparse%dim1),sparse%ia(2:sparse%dim1+1),&
@@ -393,6 +408,12 @@ module subroutine multgenm_csr(sparse,alpha,trans,x,val,y)
                  x,size(x,1),&
                  val,y,size(y,1))
 #endif
+#else
+ call csrmm(trans,sparse%dim1,size(y,2),sparse%dim2,&
+            alpha,matdescra,sparse%a,sparse%ja,sparse%ia(1:sparse%dim1),sparse%ia(2:sparse%dim1+1),&
+            x,size(x,1),&
+            val,y,size(y,1))
+#endif
 
 end subroutine
 
@@ -1193,22 +1214,30 @@ module subroutine isolve_crs(sparse,x,y)
  !$ t2=omp_get_wtime()
 #endif
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrtrsv('U','T','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
 #else
  call mkl_dcsrtrsv('U','T','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
 #endif
+#else
+  call csrtrsv('U','T','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
+#endif
 
 #if (_VERBOSE>0)
  !$ write(sparse%unlog,'(1x,a,t30,a,f0.5)')'ISOLVE CRS 1st triangular solve',': Elapsed time (s) = ',omp_get_wtime()-t2
  !$ t2=omp_get_wtime()
 #endif
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrtrsv('U','N','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
 #else
  call mkl_dcsrtrsv('U','N','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
 #endif
+#else
+ call csrtrsv('U','N','N',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
+#endif
 
 #if (_VERBOSE>0)
  !$ write(sparse%unlog,'(1x,a,t30,a,f0.5)')'ISOLVE CRS 2nd triangular solve',': Elapsed time (s) = ',omp_get_wtime()-t2
@@ -1253,11 +1282,15 @@ module subroutine solveldlt_s_crs(sparse,x,y)
  !$ t2=omp_get_wtime()
 #endif
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrtrsv('U','T','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
 #else
  call mkl_dcsrtrsv('U','T','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
 #endif
+#else
+ call csrtrsv('U','T','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x_,x)
+#endif
 #if (_VERBOSE>0)
  !$ write(sparse%unlog,'(1x,a,t30,a,f0.5)')'SOLVE LDLt CRS 1st triangular solve',': Elapsed time (s) = ',omp_get_wtime()-t2
  !$ t2=omp_get_wtime()
@@ -1275,11 +1308,15 @@ module subroutine solveldlt_s_crs(sparse,x,y)
  !$ t2=omp_get_wtime()
 #endif
 
+#if (_MKL_SPARSE==1)
 #if(_DP==0)
  call mkl_scsrtrsv('U','N','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
 #else
  call mkl_dcsrtrsv('U','N','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
 #endif
+#else
+ call csrtrsv('U','N','U',sparse%getdim(1),sparse%a,sparse%ia,sparse%ja,x,x_)
+#endif
 
 #if (_VERBOSE>0)
  !$ write(sparse%unlog,'(1x,a,t30,a,f0.5)')'SOLVE LDLt CRS 2nd triangular solve',': Elapsed time (s) = ',omp_get_wtime()-t2