diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index eacd2b7..b4e90ee 100644
--- a/.github/workflows/CI.yml
+++ b/.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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0e0c59b..67c9318 100644
--- a/CMakeLists.txt
+++ b/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}"
diff --git a/Makefile b/Makefile
index ded44ef..98e6fc5 100644
--- a/Makefile
+++ b/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 =
 
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 737fa4e..1d257cf 100644
--- a/src/CMakeLists.txt
+++ b/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"
diff --git a/src/Makefile b/src/Makefile
index b18d131..68fbb9e 100644
--- a/src/Makefile
+++ b/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
diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index c9e32ec..c182039 100644
--- a/src/modsparse_crs.f90
+++ b/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
diff --git a/src/modsparse_helpers.f90 b/src/modsparse_helpers.f90
new file mode 100644
index 0000000..9526a6e
--- /dev/null
+++ b/src/modsparse_helpers.f90
@@ -0,0 +1,238 @@
+!> Module containing various helpers
+!> Most likely all of them inefficient but should do the job
+
+module modsparse_helpers
+#if (_DP==0)
+ use iso_fortran_env,only:output_unit,int32,int64,real32,real64,wp=>real32
+#else
+ use iso_fortran_env,only:output_unit,int32,int64,real32,real64,wp=>real64
+#endif
+ !$ use omp_lib
+ implicit none(type, external)
+ private
+ public :: csrsymv
+ public :: csrmm
+ public :: csrmv
+ public :: csrtrsv
+
+ contains
+
+! Symmetric CSR matrix-vector product: y = A*x
+subroutine csrsymv(uplo, n, a, ia, ja, x, y)
+ character(len=1), intent(in) :: uplo
+ integer, intent(in) :: n
+ integer(kind=int32), intent(in) :: ia(:)
+ integer(kind=int32), intent(in) :: ja(:)
+ real(kind=wp), intent(in) :: a(:)
+ real(kind=wp), intent(in) :: x(:)
+ real(kind=wp), intent(out) :: y(:)
+
+ integer :: i, j
+
+ y = 0._wp
+
+ select case(uplo)
+ case ('U', 'u')
+  do i = 1, n
+   do j = ia(i), ia(i+1)-1
+    y(i) = y(i) + a(j) * x(ja(j))
+    if (ja(j) /= i) y(ja(j)) = y(ja(j)) + a(j) * x(i)
+   enddo
+  enddo
+ case ('L', 'l')
+  do i = 1, n
+   do j = ia(i), ia(i+1)-1
+    y(i) = y(i) + a(j) * x(ja(j))
+    if (ja(j) /= i) y(ja(j)) = y(ja(j)) + a(j) * x(i)
+   enddo
+  enddo
+ case default
+  error stop 'csrsymv: unsupported uplo'
+ end select
+
+end subroutine csrsymv
+
+! CSR matrix-dense matrix product: c = beta*c + alpha*op(A)*b
+subroutine csrmm(transa, m, n, k, alpha, matdescra, a, ja, pntrb, pntre, b, ldb, beta, c, ldc)
+ character(len=1), intent(in) :: transa
+ integer, intent(in) :: m
+ integer, intent(in) :: n
+ integer, intent(in) :: k
+ real(kind=wp), intent(in) :: alpha
+ character(len=1), intent(in) :: matdescra(6)
+ real(kind=wp), intent(in) :: a(:)
+ integer(kind=int32), intent(in) :: ja(:)
+ integer(kind=int32), intent(in) :: pntrb(:)
+ integer(kind=int32), intent(in) :: pntre(:)
+ real(kind=wp), intent(in) :: b(:,:)
+ integer, intent(in) :: ldb
+ real(kind=wp), intent(in) :: beta
+ real(kind=wp), intent(inout) :: c(:,:)
+ integer, intent(in) :: ldc
+
+ integer :: i, j, p
+
+ if (transa == 'N' .or. transa == 'n') then
+
+  c(1:m, 1:n) = beta * c(1:m, 1:n)
+
+  select case (matdescra(1))
+  case ('G', 'g', 'T', 't')
+   do j = 1, n
+    do i = 1, m
+     do p = pntrb(i), pntre(i)-1
+      c(i,j) = c(i,j) + alpha * a(p) * b(ja(p),j)
+     end do
+    end do
+   end do
+  case ('S', 's')
+   do j = 1, n
+    do i = 1, m
+     do p = pntrb(i), pntre(i)-1
+      c(i,j) = c(i,j) + alpha * a(p) * b(ja(p),j)
+      if (ja(p) /= i) c(ja(p),j) = c(ja(p),j) + alpha * a(p) * b(i,j)
+     end do
+    end do
+   end do
+  case default
+   error stop 'csrmm: unsupported matdescra(1)'
+  end select
+
+ else
+
+  select case (matdescra(1))
+  case ('S', 's')
+   ! Symmetric: A^T = A
+   c(1:m, 1:n) = beta * c(1:m, 1:n)
+   do j = 1, n
+    do i = 1, m
+     do p = pntrb(i), pntre(i)-1
+      c(i,j) = c(i,j) + alpha * a(p) * b(ja(p),j)
+      if (ja(p) /= i) c(ja(p),j) = c(ja(p),j) + alpha * a(p) * b(i,j)
+     end do
+    end do
+   end do
+  case ('G', 'g', 'T', 't')
+   c(1:k, 1:n) = beta * c(1:k, 1:n)
+   do j = 1, n
+    do i = 1, m
+     do p = pntrb(i), pntre(i)-1
+      c(ja(p),j) = c(ja(p),j) + alpha * a(p) * b(i,j)
+     end do
+    end do
+   end do
+  case default
+   error stop 'csrmm: unsupported matdescra(1)'
+  end select
+
+ end if
+
+end subroutine csrmm
+
+! General CSR matrix-vector product: y = beta*y + alpha*op(A)*x
+subroutine csrmv(transa, m, k, alpha, matdescra, val, indx, pntrb, pntre, x, beta, y)
+ character(len=1), intent(in) :: transa
+ integer, intent(in) :: m
+ integer, intent(in) :: k
+ real(kind=wp), intent(in) :: alpha
+ real(kind=wp), intent(in) :: beta
+ character(len=1), intent(in) :: matdescra(:)
+ real(kind=wp), intent(in) :: val(:)
+ real(kind=wp), intent(in) :: x(:)
+ integer(kind=int32), intent(in) :: indx(:)
+ integer(kind=int32), intent(in) :: pntrb(:)
+ integer(kind=int32), intent(in) :: pntre(:)
+ real(kind=wp), intent(inout) :: y(:)
+
+ integer :: i, j
+
+ if (transa == 'N' .or. transa == 'n') then
+  y(1:m) = beta * y(1:m)
+  select case (matdescra(1))
+  case ('G', 'g', 'T', 't')
+   do i = 1, m
+    do j = pntrb(i), pntre(i)-1
+     y(i) = y(i) + alpha * val(j) * x(indx(j))
+    end do
+   end do
+  case ('S', 's')
+   do i = 1, m
+    do j = pntrb(i), pntre(i)-1
+     y(i) = y(i) + alpha * val(j) * x(indx(j))
+     if (indx(j) /= i) y(indx(j)) = y(indx(j)) + alpha * val(j) * x(i)
+    end do
+   end do
+  case default
+   error stop 'csrmv: unsupported matdescra(1)'
+  end select
+ else
+  select case (matdescra(1))
+  case ('G', 'g', 'T', 't')
+   y(1:k) = beta * y(1:k)
+   do i = 1, m
+    do j = pntrb(i), pntre(i)-1
+     y(indx(j)) = y(indx(j)) + alpha * val(j) * x(i)
+    end do
+   end do
+  case ('S', 's')
+   y(1:m) = beta * y(1:m)
+   do i = 1, m
+    do j = pntrb(i), pntre(i)-1
+     y(i) = y(i) + alpha * val(j) * x(indx(j))
+     if (indx(j) /= i) y(indx(j)) = y(indx(j)) + alpha * val(j) * x(i)
+    end do
+   end do
+  case default
+   error stop 'csrmv: unsupported matdescra(1)'
+  end select
+ end if
+
+end subroutine csrmv
+
+! Upper triangular CSR solve: op(A)*y = x
+! Diagonal is first entry in each row (sorted CSR assumed)
+subroutine csrtrsv(uplo, transa, diag, m, a, ia, ja, x, y)
+ character(len=1), intent(in) :: uplo
+ character(len=1), intent(in) :: transa
+ character(len=1), intent(in) :: diag
+ integer, intent(in) :: m
+ real(kind=wp), intent(in) :: a(:)
+ integer(kind=int32), intent(in) :: ia(:)
+ integer(kind=int32), intent(in) :: ja(:)
+ real(kind=wp), intent(in) :: x(:)
+ real(kind=wp), intent(out) :: y(:)
+
+ integer :: i, k
+ logical :: ldiag
+
+ y = x
+
+ select case(uplo)
+ case ('U', 'u')
+  ldiag = (diag == 'N' .or. diag == 'n')
+
+  if (transa == 'N' .or. transa == 'n') then
+   ! Solve U*y = x; backward substitution
+   do i = m, 1, -1
+    do k = ia(i)+1, ia(i+1)-1
+     y(i) = y(i) - a(k) * y(ja(k))
+    end do
+    if (ldiag) y(i) = y(i) / a(ia(i))
+   end do
+  else
+   ! Solve U^T*y = x; forward substitution
+   do i = 1, m
+    if (ldiag) y(i) = y(i) / a(ia(i))
+    do k = ia(i)+1, ia(i+1)-1
+     y(ja(k)) = y(ja(k)) - a(k) * y(i)
+    end do
+   end do
+  end if
+
+ case default
+  error stop 'internal error csrtrsv'
+ end select
+
+end subroutine csrtrsv
+
+end module modsparse_helpers