From b4500ca258fa514b2b1e63553eb62b7bdf80be3c Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 13:29:56 +0100
Subject: [PATCH 1/8] modify sort_crs for later

---
 src/modsparse_crs.f90 | 37 +++++++++++++++++++++++++++----------
 1 file changed, 27 insertions(+), 10 deletions(-)

diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index ccc1290..b88bb02 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -1163,13 +1163,34 @@ module subroutine sort_crs(sparse)
   return
  endif
 
- do k=1,sparse%dim1
-  n=sparse%ia(k+1)-sparse%ia(k)
+ call sort_crs_ija(sparse%dim1, sparse%ia, sparse%ja, sparse%a)
+
+ call sparse%setsorted(.true.)
+
+end subroutine
+
+module subroutine sort_crs_ija(dim1, ia, ja, a)
+ ! sort vectors ja and a by increasing order
+ integer(kind=int32), intent(in) :: dim1
+ integer(kind=int32), intent(in) :: ia(:)
+ integer(kind=int32), intent(inout) :: ja(:)
+ real(kind=wp), intent(inout) :: a(:)
+
+ integer(kind=int32)::endd,i,j,k,n,start,stkpnt
+ integer(kind=int32)::d1,d2,d3,dmnmx,tmp
+ integer(kind=int32)::stack(2,32)
+ integer(kind=int32),allocatable::d(:)
+ integer(kind=int32),parameter::select=20
+ real(kind=wp)::umnmx,tmpu
+ real(kind=wp),allocatable::u(:)
+
+ do k = 1, dim1
+  n = ia(k+1) - ia(k)
   if(n.gt.1)then
    allocate(d(n),u(n))
    !copy of the vector to be sorted
-   d=sparse%ja(sparse%ia(k):sparse%ia(k+1)-1)
-   u=sparse%a(sparse%ia(k):sparse%ia(k+1)-1)
+   d=ja(ia(k):ia(k+1)-1)
+   u=a(ia(k):ia(k+1)-1)
    !sort the vectors
    !from dlasrt.f !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !Quick return if possible
@@ -1255,17 +1276,13 @@ module subroutine sort_crs(sparse)
    IF( stkpnt > 0 ) GO TO 10
    !end from dlasrt.f !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !copy back the sorted vectors
-   sparse%ja(sparse%ia(k):sparse%ia(k+1)-1)=d
-   sparse%a(sparse%ia(k):sparse%ia(k+1)-1)=u
+   ja(ia(k):ia(k+1)-1)=d
+   a(ia(k):ia(k+1)-1)=u
    deallocate(d,u)
   endif
  enddo
 
- call sparse%setsorted(.true.)
-
 end subroutine
 
 
-
-
 end submodule

From f3bd3a138ffbfef07a202dc5bff0ead69a5b0c25 Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 15:41:00 +0100
Subject: [PATCH 2/8] add the subroutine extractcrs_fromcoo_int32

---
 src/modsparse.f90     | 128 ++++++++++++++++++++++++++++++++++++++++++
 src/modsparse_crs.f90 |   5 +-
 test/modtest_coo.f90  |  47 +++++++++++++++-
 3 files changed, 174 insertions(+), 6 deletions(-)

diff --git a/src/modsparse.f90 b/src/modsparse.f90
index a1dabfc..ad270e5 100644
--- a/src/modsparse.f90
+++ b/src/modsparse.f90
@@ -23,6 +23,7 @@ module modsparse
  public::gen_sparse
  public::llsparse,coosparse,crssparse
  public::assignment(=)
+ public::extractcrs
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!GEN!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!aaa
  integer(kind=int32),parameter::typegen=1,typecoo=10,typecrs=20,typell=30
@@ -861,6 +862,21 @@ module function getmem_metisgraph(metis) result(getmem)
                   ,convertfromcrstocoo,convertfromcrstoll
  end interface
 
+ interface
+  module subroutine sort_crs_ija_int32(dim1, ia, ja, a)
+   ! sort vectors ja and a by increasing order
+   implicit none
+   integer(kind=int32), intent(in) :: dim1
+   integer(kind=int32), intent(in) :: ia(:)
+   integer(kind=int32), intent(inout) :: ja(:)
+   real(kind=wp), intent(inout) :: a(:)
+  end subroutine
+ end interface
+
+ interface extractcrs
+  module procedure extractcrs_fromcoo_int32
+ end interface
+
 contains
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!COO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!aaa
@@ -1639,4 +1655,116 @@ subroutine convertfromcrstometisgraph(metis,sparse)
  enddo
 
 end subroutine
+
+!EXTRACT DATA FROM SPARSE ARRAYS
+subroutine extractcrs_fromcoo_int32(sparse, dim1, dim2, ia, ja, a, perm)
+ type(coosparse),intent(in)::sparse
+ integer(kind=int32) :: dim1, dim2
+ integer(kind=int32), allocatable :: ia(:)
+ integer(kind=int32), allocatable :: ja(:)
+ real(kind=wp), allocatable :: a(:)
+ integer(kind=int32), allocatable, optional :: perm(:)
+
+ integer(kind=int32)::i,nel,row,col
+ integer(kind=int32),allocatable::rowpos(:)
+ integer(kind=int64)::i8
+ logical::lsquare
+
+ if(sparse%nonzero().ge.huge(ia)/2)then
+  write(sparse%unlog,'(a)')' ERROR: impossible conversion due to a too large number of non-zero elements'
+  error stop
+ endif
+
+ !Condition: all rows contain at least one element (diagonal element if square or one dummy entry in the last column if needed)
+
+ !Number of elements=number of rows+number off-diagonal elements
+ !from sparse
+ lsquare=sparse%issquare()
+
+ allocate(rowpos(sparse%dim1))
+ rowpos=0
+
+ if(lsquare)then
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.ne.0.and.row.ne.sparse%ij(2,i8))then
+    rowpos(row)=rowpos(row)+1
+   endif
+  enddo
+ else
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.ne.0.and.sparse%ij(2,i8).ne.sparse%dim2)then
+    rowpos(row)=rowpos(row)+1
+   endif
+  enddo
+ endif
+
+ nel=sparse%dim1+sum(rowpos)
+
+ dim1 = sparse%dim1
+ dim2 = sparse%dim2
+ allocate(ia(sparse%dim1+1), source = 0)
+ ia(sparse%dim1+1) = -nel
+ allocate(ja(nel), source = 0)
+ allocate(a(nel), source = 0._wp)
+
+ if(allocated(sparse%perm))allocate(perm, source = sparse%perm)
+
+ !determine the number of non-zero off-diagonal elements per row
+ ia(2:sparse%dim1+1) = rowpos
+
+ !accumulate the number of elements and add diagonal elements
+ ia(1) = 1
+ if(lsquare)then
+  do i=1,sparse%dim1
+   ia(i+1) = ia(i+1) + 1 + ia(i)
+   ja(ia(i)) = i   !set diagonal element for the case it would not be present
+  enddo
+ else
+  do i=1,sparse%dim1
+   ia(i+1) = ia(i+1) + 1 + ia(i)
+   ja(ia(i)) = sparse%dim2   !set last element for the case it would not be present
+  enddo
+ endif
+
+ !add the non-zero elements to crs
+ !allocate(rowpos(dim1))
+ rowpos = ia(1:sparse%dim1)
+ if(lsquare)then
+  do i8 = 1_int64, sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.gt.0)then
+    col=sparse%ij(2,i8)
+    if(col.eq.row)then
+     a(ia(row)) = sparse%a(i8)
+    else
+     rowpos(row) = rowpos(row)+1
+     ja(rowpos(row)) = col
+     a(rowpos(row)) = sparse%a(i8)
+    endif
+   endif
+  enddo
+ else
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.gt.0)then
+    col=sparse%ij(2,i8)
+    if(col.eq.sparse%dim2)then
+     a(ia(row)) = sparse%a(i8)
+    else
+     rowpos(row) = rowpos(row)+1
+     ja(rowpos(row)) = col
+     a(rowpos(row)) = sparse%a(i8)
+    endif
+   endif
+  enddo
+ endif
+
+ deallocate(rowpos)
+
+ call sort_crs_ija_int32(dim1, ia, ja, a)
+
+end subroutine
+
 end module
diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index b88bb02..2bbb51e 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -1163,13 +1163,13 @@ module subroutine sort_crs(sparse)
   return
  endif
 
- call sort_crs_ija(sparse%dim1, sparse%ia, sparse%ja, sparse%a)
+ call sort_crs_ija_int32(sparse%dim1, sparse%ia, sparse%ja, sparse%a)
 
  call sparse%setsorted(.true.)
 
 end subroutine
 
-module subroutine sort_crs_ija(dim1, ia, ja, a)
+subroutine sort_crs_ija_int32(dim1, ia, ja, a)
  ! sort vectors ja and a by increasing order
  integer(kind=int32), intent(in) :: dim1
  integer(kind=int32), intent(in) :: ia(:)
@@ -1284,5 +1284,4 @@ module subroutine sort_crs_ija(dim1, ia, ja, a)
 
 end subroutine
 
-
 end submodule
diff --git a/test/modtest_coo.f90 b/test/modtest_coo.f90
index 305a51c..8776fff 100644
--- a/test/modtest_coo.f90
+++ b/test/modtest_coo.f90
@@ -1,11 +1,11 @@
 module modtest_coo
 #if (_DP==0)
- use, intrinsic :: iso_fortran_env, only: int64, output_unit, wp => real32
+ use, intrinsic :: iso_fortran_env, only: int32, int64, output_unit, wp => real32
 #else
- use, intrinsic :: iso_fortran_env, only: int64, output_unit, wp => real64
+ use, intrinsic :: iso_fortran_env, only: int32, int64, output_unit, wp => real64
 #endif
  use testdrive, only: new_unittest, unittest_type, error_type, check
- use modsparse, only: coosparse
+ use modsparse, only: coosparse, extractcrs
  use modtest_common, only: tol_wp, verbose, ia, ja, a, aspsd&
                        , addval => addval_coo, getmat, matcheck, printmat
  implicit none
@@ -35,6 +35,7 @@ subroutine collect_coo(testsuite)
     , new_unittest("coo diag vect lupper", test_diag_vect_lupper) &
     , new_unittest("coo ncol diag vect", test_ncol_diag_vect) &
     , new_unittest("coo ncol diag vect lupper", test_ncol_diag_vect_lupper) &
+    , new_unittest("coo extract lupper_sym", test_extract_lupper_sym) &
     , new_unittest("coo get", test_get) &
     , new_unittest("coo get nel", test_get_nel) &
     , new_unittest("coo get lupper", test_get_lupper) &
@@ -464,6 +465,46 @@ subroutine test_ncol_diag_vect_lupper(error)
 
 end subroutine
 
+!EXTRACT
+subroutine test_extract_lupper_sym(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nrow = 5
+ integer, parameter :: ncol = nrow
+ logical, parameter :: lvalid(size(ia)) = ia.le.nrow .and. ja.le.ncol .and. ia .le. ja
+
+ type(coosparse) :: coo
+ integer(kind=int32) :: dim1, dim2
+ integer(kind=int32), allocatable :: iatmp(:), jatmp(:)
+ real(kind=wp), allocatable :: atmp(:)
+
+ coo = coosparse(nrow, lupper = .true., unlog = sparse_unit)
+ call coo%setsymmetric()
+
+ call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, a)
+
+ call extractcrs(coo, dim1, dim2, iatmp, jatmp, atmp)
+
+ call check(error, dim1 == nrow, 'extract lupper_sym, dim1')
+ if(allocated(error))return
+
+ call check(error, dim2 == ncol, 'extract lupper_sym, dim2')
+ if(allocated(error))return
+
+ call check(error, all(iatmp == [1, 5, 7, 9, 10 , 11]),  'extract lupper_sym, ia')
+ if(allocated(error))return
+
+ call check(error, all(jatmp == [1, 3, 4, 5, 2, 5, 3, 4, 4, 5])&
+             ,  'extract lupper_sym, ja')
+ if(allocated(error))return
+
+ call check(error, all(atmp == [11._wp, 13._wp, 14._wp, 15._wp, 0._wp&
+             , 25._wp, 33._wp, 34._wp, 44._wp, 55._wp]) &
+             ,  'extract lupper_sym, a')
+ if(allocated(error))return
+
+end subroutine
+
 !GET
 subroutine test_get(error)
  type(error_type), allocatable, intent(out) :: error

From 6a4f35725eb96d0aab55a6eddd752750dcd0c562 Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 15:54:14 +0100
Subject: [PATCH 3/8] extractcrs_fromcoo_int64

---
 src/modsparse.f90     | 119 ++++++++++++++++++++++++++++++++++++++++++
 src/modsparse_crs.f90 | 114 ++++++++++++++++++++++++++++++++++++++++
 test/modtest_coo.f90  |  40 ++++++++++++++
 3 files changed, 273 insertions(+)

diff --git a/src/modsparse.f90 b/src/modsparse.f90
index ad270e5..61a89a1 100644
--- a/src/modsparse.f90
+++ b/src/modsparse.f90
@@ -871,10 +871,19 @@ module subroutine sort_crs_ija_int32(dim1, ia, ja, a)
    integer(kind=int32), intent(inout) :: ja(:)
    real(kind=wp), intent(inout) :: a(:)
   end subroutine
+  module subroutine sort_crs_ija_int64(dim1, ia, ja, a)
+   ! sort vectors ja and a by increasing order
+   implicit none
+   integer(kind=int64), intent(in) :: dim1
+   integer(kind=int64), intent(in) :: ia(:)
+   integer(kind=int64), intent(inout) :: ja(:)
+   real(kind=wp), intent(inout) :: a(:)
+  end subroutine
  end interface
 
  interface extractcrs
   module procedure extractcrs_fromcoo_int32
+  module procedure extractcrs_fromcoo_int64
  end interface
 
 contains
@@ -1767,4 +1776,114 @@ subroutine extractcrs_fromcoo_int32(sparse, dim1, dim2, ia, ja, a, perm)
 
 end subroutine
 
+subroutine extractcrs_fromcoo_int64(sparse, dim1, dim2, ia, ja, a, perm)
+ type(coosparse),intent(in)::sparse
+ integer(kind=int64) :: dim1, dim2
+ integer(kind=int64), allocatable :: ia(:)
+ integer(kind=int64), allocatable :: ja(:)
+ real(kind=wp), allocatable :: a(:)
+ integer(kind=int64), allocatable, optional :: perm(:)
+
+ integer(kind=int32)::i,row,col
+ integer(kind=int64)::nel
+ integer(kind=int64),allocatable::rowpos(:)
+ integer(kind=int64)::i8
+ logical::lsquare
+
+ if(sparse%nonzero().ge.huge(ia)/2)then
+  write(sparse%unlog,'(a)')' ERROR: impossible conversion due to a too large number of non-zero elements'
+  error stop
+ endif
+
+ !Condition: all rows contain at least one element (diagonal element if square or one dummy entry in the last column if needed)
+
+ !Number of elements=number of rows+number off-diagonal elements
+ !from sparse
+ lsquare=sparse%issquare()
+
+ allocate(rowpos(sparse%dim1))
+ rowpos=0
+
+ if(lsquare)then
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.ne.0.and.row.ne.sparse%ij(2,i8))then
+    rowpos(row)=rowpos(row)+1
+   endif
+  enddo
+ else
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.ne.0.and.sparse%ij(2,i8).ne.sparse%dim2)then
+    rowpos(row)=rowpos(row)+1
+   endif
+  enddo
+ endif
+
+ nel=sparse%dim1+sum(rowpos)
+
+ dim1 = sparse%dim1
+ dim2 = sparse%dim2
+ allocate(ia(sparse%dim1+1), source = 0_int64)
+ ia(sparse%dim1+1) = -nel
+ allocate(ja(nel), source = 0_int64)
+ allocate(a(nel), source = 0._wp)
+
+ if(allocated(sparse%perm))allocate(perm, source = int(sparse%perm, int64))
+
+ !determine the number of non-zero off-diagonal elements per row
+ ia(2:sparse%dim1+1) = rowpos
+
+ !accumulate the number of elements and add diagonal elements
+ ia(1) = 1
+ if(lsquare)then
+  do i=1,sparse%dim1
+   ia(i+1) = ia(i+1) + 1 + ia(i)
+   ja(ia(i)) = i   !set diagonal element for the case it would not be present
+  enddo
+ else
+  do i=1,sparse%dim1
+   ia(i+1) = ia(i+1) + 1 + ia(i)
+   ja(ia(i)) = sparse%dim2   !set last element for the case it would not be present
+  enddo
+ endif
+
+ !add the non-zero elements to crs
+ !allocate(rowpos(dim1))
+ rowpos = ia(1:sparse%dim1)
+ if(lsquare)then
+  do i8 = 1_int64, sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.gt.0)then
+    col=sparse%ij(2,i8)
+    if(col.eq.row)then
+     a(ia(row)) = sparse%a(i8)
+    else
+     rowpos(row) = rowpos(row)+1
+     ja(rowpos(row)) = col
+     a(rowpos(row)) = sparse%a(i8)
+    endif
+   endif
+  enddo
+ else
+  do i8=1_int64,sparse%nel
+   row=sparse%ij(1,i8)
+   if(row.gt.0)then
+    col=sparse%ij(2,i8)
+    if(col.eq.sparse%dim2)then
+     a(ia(row)) = sparse%a(i8)
+    else
+     rowpos(row) = rowpos(row)+1
+     ja(rowpos(row)) = col
+     a(rowpos(row)) = sparse%a(i8)
+    endif
+   endif
+  enddo
+ endif
+
+ deallocate(rowpos)
+
+ call sort_crs_ija_int64(dim1, ia, ja, a)
+
+end subroutine
 end module
diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index 2bbb51e..7420e88 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -1284,4 +1284,118 @@ subroutine sort_crs_ija_int32(dim1, ia, ja, a)
 
 end subroutine
 
+subroutine sort_crs_ija_int64(dim1, ia, ja, a)
+ ! sort vectors ja and a by increasing order
+ integer(kind=int64), intent(in) :: dim1
+ integer(kind=int64), intent(in) :: ia(:)
+ integer(kind=int64), intent(inout) :: ja(:)
+ real(kind=wp), intent(inout) :: a(:)
+
+ integer(kind=int64)::endd,i,j,k,n,start,stkpnt
+ integer(kind=int64)::d1,d2,d3,dmnmx,tmp
+ integer(kind=int64)::stack(2,32)
+ integer(kind=int64),allocatable::d(:)
+ integer(kind=int64),parameter::select=20
+ real(kind=wp)::umnmx,tmpu
+ real(kind=wp),allocatable::u(:)
+
+ do k = 1, dim1
+  n = ia(k+1) - ia(k)
+  if(n.gt.1)then
+   allocate(d(n),u(n))
+   !copy of the vector to be sorted
+   d=ja(ia(k):ia(k+1)-1)
+   u=a(ia(k):ia(k+1)-1)
+   !sort the vectors
+   !from dlasrt.f !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   !Quick return if possible
+   stkpnt = 1
+   stack( 1, 1 ) = 1
+   stack( 2, 1 ) = n
+   10 start = stack( 1, stkpnt )
+   endd = stack( 2, stkpnt )
+   stkpnt = stkpnt - 1
+   IF( endd-start <= select .AND. endd-start > 0 ) THEN
+   !Do Insertion sort on D( START:ENDD )
+   !Sort into increasing order
+     DO i = start + 1, endd
+      DO j = i, start + 1, -1
+       IF( d( j ) < d( j-1 ) ) THEN
+         dmnmx = d( j )
+         d( j ) = d( j-1 )
+         d( j-1 ) = dmnmx
+         umnmx = u( j )
+         u( j ) = u( j-1 )
+         u( j-1 ) = umnmx
+       ELSE
+         CYCLE
+       END IF
+      END DO
+     END DO
+   ELSE IF( endd-start > select ) THEN
+   !Partition D( START:ENDD ) and stack parts, largest one first
+   !Choose partition entry as median of 3
+    d1 = d( start )
+    d2 = d( endd )
+    i = ( start + endd ) / 2
+    d3 = d( i )
+    IF( d1 < d2 ) THEN
+      IF( d3 < d1 ) THEN
+        dmnmx = d1
+      ELSE IF( d3 < d2 ) THEN
+        dmnmx = d3
+      ELSE
+        dmnmx = d2
+      END IF
+    ELSE
+      IF( d3 < d2 ) THEN
+        dmnmx = d2
+      ELSE IF( d3 < d1 ) THEN
+        dmnmx = d3
+      ELSE
+        dmnmx = d1
+      END IF
+    END IF
+    !Sort into increasing order
+     i = start - 1
+     j = endd + 1
+     90 j = j - 1
+     IF( d( j ) > dmnmx ) GO TO 90
+     110 i = i + 1
+     IF( d( i ) < dmnmx ) GO TO 110
+     IF( i < j ) THEN
+       tmp = d( i )
+       d( i ) = d( j )
+       d( j ) = tmp
+       tmpu = u( i )
+       u( i ) = u( j )
+       u( j ) = tmpu
+       GO TO 90
+     END IF
+     IF( j-start > endd-j-1 ) THEN
+       stkpnt = stkpnt + 1
+       stack( 1, stkpnt ) = start
+       stack( 2, stkpnt ) = j
+       stkpnt = stkpnt + 1
+       stack( 1, stkpnt ) = j + 1
+       stack( 2, stkpnt ) = endd
+     ELSE
+       stkpnt = stkpnt + 1
+       stack( 1, stkpnt ) = j + 1
+       stack( 2, stkpnt ) = endd
+       stkpnt = stkpnt + 1
+       stack( 1, stkpnt ) = start
+       stack( 2, stkpnt ) = j
+     END IF
+   END IF
+   IF( stkpnt > 0 ) GO TO 10
+   !end from dlasrt.f !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   !copy back the sorted vectors
+   ja(ia(k):ia(k+1)-1)=d
+   a(ia(k):ia(k+1)-1)=u
+   deallocate(d,u)
+  endif
+ enddo
+
+end subroutine
 end submodule
diff --git a/test/modtest_coo.f90 b/test/modtest_coo.f90
index 8776fff..06b8f6f 100644
--- a/test/modtest_coo.f90
+++ b/test/modtest_coo.f90
@@ -36,6 +36,7 @@ subroutine collect_coo(testsuite)
     , new_unittest("coo ncol diag vect", test_ncol_diag_vect) &
     , new_unittest("coo ncol diag vect lupper", test_ncol_diag_vect_lupper) &
     , new_unittest("coo extract lupper_sym", test_extract_lupper_sym) &
+    , new_unittest("coo extract lupper_sym_int64", test_extract_lupper_sym_int64) &
     , new_unittest("coo get", test_get) &
     , new_unittest("coo get nel", test_get_nel) &
     , new_unittest("coo get lupper", test_get_lupper) &
@@ -505,6 +506,45 @@ subroutine test_extract_lupper_sym(error)
 
 end subroutine
 
+subroutine test_extract_lupper_sym_int64(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nrow = 5
+ integer, parameter :: ncol = nrow
+ logical, parameter :: lvalid(size(ia)) = ia.le.nrow .and. ja.le.ncol .and. ia .le. ja
+
+ type(coosparse) :: coo
+ integer(kind=int64) :: dim1, dim2
+ integer(kind=int64), allocatable :: iatmp(:), jatmp(:)
+ real(kind=wp), allocatable :: atmp(:)
+
+ coo = coosparse(nrow, lupper = .true., unlog = sparse_unit)
+ call coo%setsymmetric()
+
+ call addval(coo, coo%getdim(1), coo%getdim(2), ia, ja, a)
+
+ call extractcrs(coo, dim1, dim2, iatmp, jatmp, atmp)
+
+ call check(error, dim1 == nrow, 'extract lupper_sym_int64, dim1')
+ if(allocated(error))return
+
+ call check(error, dim2 == ncol, 'extract lupper_sym_int64, dim2')
+ if(allocated(error))return
+
+ call check(error, all(iatmp == [integer(int64) :: 1, 5, 7, 9, 10 , 11]),  'extract lupper_sym_int64, ia')
+ if(allocated(error))return
+
+ call check(error, all(jatmp == [integer(int64) :: 1, 3, 4, 5, 2, 5, 3, 4, 4, 5])&
+             ,  'extract lupper_sym_int64, ja')
+ if(allocated(error))return
+
+ call check(error, all(atmp == [11._wp, 13._wp, 14._wp, 15._wp, 0._wp&
+             , 25._wp, 33._wp, 34._wp, 44._wp, 55._wp]) &
+             ,  'extract lupper_sym_int64, a')
+ if(allocated(error))return
+
+end subroutine
+
 !GET
 subroutine test_get(error)
  type(error_type), allocatable, intent(out) :: error

From b5cd36281afab4953ead35bbe098d262cf0c25e5 Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 16:55:29 +0100
Subject: [PATCH 4/8] fix issue with declaring sort_crs_int64

---
 src/modsparse_crs.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index 7420e88..dcf16be 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -1169,7 +1169,7 @@ module subroutine sort_crs(sparse)
 
 end subroutine
 
-subroutine sort_crs_ija_int32(dim1, ia, ja, a)
+module subroutine sort_crs_ija_int32(dim1, ia, ja, a)
  ! sort vectors ja and a by increasing order
  integer(kind=int32), intent(in) :: dim1
  integer(kind=int32), intent(in) :: ia(:)
@@ -1284,7 +1284,7 @@ subroutine sort_crs_ija_int32(dim1, ia, ja, a)
 
 end subroutine
 
-subroutine sort_crs_ija_int64(dim1, ia, ja, a)
+module subroutine sort_crs_ija_int64(dim1, ia, ja, a)
  ! sort vectors ja and a by increasing order
  integer(kind=int64), intent(in) :: dim1
  integer(kind=int64), intent(in) :: ia(:)

From 864796a2aa12bfff2c1eb02d58d595efd2101f75 Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 17:57:42 +0100
Subject: [PATCH 5/8] fix possible bug

---
 src/modsparse.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/modsparse.f90 b/src/modsparse.f90
index 61a89a1..9ce7e3c 100644
--- a/src/modsparse.f90
+++ b/src/modsparse.f90
@@ -1718,7 +1718,7 @@ subroutine extractcrs_fromcoo_int32(sparse, dim1, dim2, ia, ja, a, perm)
  allocate(ja(nel), source = 0)
  allocate(a(nel), source = 0._wp)
 
- if(allocated(sparse%perm))allocate(perm, source = sparse%perm)
+ if(allocated(sparse%perm).and.present(perm))allocate(perm, source = int(sparse%perm, int64))
 
  !determine the number of non-zero off-diagonal elements per row
  ia(2:sparse%dim1+1) = rowpos
@@ -1829,7 +1829,7 @@ subroutine extractcrs_fromcoo_int64(sparse, dim1, dim2, ia, ja, a, perm)
  allocate(ja(nel), source = 0_int64)
  allocate(a(nel), source = 0._wp)
 
- if(allocated(sparse%perm))allocate(perm, source = int(sparse%perm, int64))
+ if(allocated(sparse%perm).and.present(perm))allocate(perm, source = int(sparse%perm, int64))
 
  !determine the number of non-zero off-diagonal elements per row
  ia(2:sparse%dim1+1) = rowpos

From f45d9dc8004cc9a8580db1a285a6bebadf5e1bc9 Mon Sep 17 00:00:00 2001
From: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
Date: Wed, 12 Jan 2022 18:00:15 +0100
Subject: [PATCH 6/8] fix issue

---
 src/modsparse.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/modsparse.f90 b/src/modsparse.f90
index 9ce7e3c..2e18c66 100644
--- a/src/modsparse.f90
+++ b/src/modsparse.f90
@@ -1718,7 +1718,7 @@ subroutine extractcrs_fromcoo_int32(sparse, dim1, dim2, ia, ja, a, perm)
  allocate(ja(nel), source = 0)
  allocate(a(nel), source = 0._wp)
 
- if(allocated(sparse%perm).and.present(perm))allocate(perm, source = int(sparse%perm, int64))
+ if(allocated(sparse%perm).and.present(perm))allocate(perm, source = sparse%perm)
 
  !determine the number of non-zero off-diagonal elements per row
  ia(2:sparse%dim1+1) = rowpos

From 54d87c563b6a83ca2d3724190fca49b308e2a26a Mon Sep 17 00:00:00 2001
From: "Vandenplas, Jeremie" <jeremie.vandenplas@wur.nl>
Date: Thu, 13 Jan 2022 09:23:54 +0100
Subject: [PATCH 7/8] unrelated change

---
 src/modsparse_crs.f90 | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index dcf16be..e150593 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -812,9 +812,12 @@ module subroutine solve_crs_vector(sparse,x,y)
 
   !initialize iparm
   call pardisoinit(sparse%pardisovar%pt,sparse%pardisovar%mtype,sparse%pardisovar%iparm)
+!  block
+!  integer::i
 !  do i=1,64
 !   write(sparse%unlog,*)'iparm',i,sparse%pardisovar%iparm(i)
 !  enddo
+!  end block
 
   !Ordering and factorization
   sparse%pardisovar%phase=12

From 1c41e5cc02a1a4b31b69fed77dad0436512f6d0d Mon Sep 17 00:00:00 2001
From: "Vandenplas, Jeremie" <jeremie.vandenplas@wur.nl>
Date: Thu, 13 Jan 2022 09:27:51 +0100
Subject: [PATCH 8/8] clean

---
 src/modsparse_crs.f90 | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/modsparse_crs.f90 b/src/modsparse_crs.f90
index e150593..23f9986 100644
--- a/src/modsparse_crs.f90
+++ b/src/modsparse_crs.f90
@@ -1154,14 +1154,6 @@ module subroutine sort_crs(sparse)
  ! sort vectors ja and a by increasing order
  class(crssparse),intent(inout)::sparse
 
- integer(kind=int32)::endd,i,j,k,n,start,stkpnt
- integer(kind=int32)::d1,d2,d3,dmnmx,tmp
- integer(kind=int32)::stack(2,32)
- integer(kind=int32),allocatable::d(:)
- integer(kind=int32),parameter::select=20
- real(kind=wp)::umnmx,tmpu
- real(kind=wp),allocatable::u(:)
-
  if(sparse%issorted())then
   return
  endif