conversion from COO format to CSR

Likun_T_ · ‎04-01-2015

Dear all,

I am using Intel MKL subroutine mkl_zcoocsr to convert a sparse matrix in COO form to CSR, but the results do not seem right. Here is my subroutine:

!ooooooooooooooooooooooooooooooooooooooooooooooooooooo

subroutine testcoocsr

      integer, dimension(:), allocatable :: irn, jcn, aval
      integer, dimension(:), allocatable :: colindex, rowindex, amtrx
      integer job(8), nnz, n, info

      nnz=6
      n=4

     allocate(irn(nnz))
     allocate(jcn(nnz))
     allocate(aval(nnz))
     allocate(colindex(nnz))
     allocate(amtrx(nnz))
     allocate(rowindex(n+1))

   job(1)=2 !from COO to CSR, column indices in CSRare sorted in increasing order within each row.
   job(2)=1 !one-based indexing for the matrix in CSR format is used.
   job(3)=1 !one-based indexing for the matrix in coordinate format is used.
   job(5)=nnz !sets number of the non-zero elements of the matrix A if job(1)=1.
   job(6)=0 !all arrays acsr, ja, ia are filled in for the output storage.

   irn(1)=1
   irn(2)=2
   irn(3)=3
   irn(4)=3
   irn(5)=4
   irn(6)=2

   jcn(1)=1
   jcn(2)=2
   jcn(3)=1
   jcn(4)=3
   jcn(5)=4
   jcn(6)=4

   aval(1)=1
   aval(2)=2
   aval(3)=3
   aval(4)=4
   aval(5)=5
   aval(6)=6

   call mkl_zcsrcoo(job, n, amtrx, colindex, rowindex, nnz, aval, irn, jcn, info)

   deallocate(irn)
     deallocate(jcn)
     deallocate(aval)
     deallocate(colindex)
     deallocate(amtrx)
     deallocate(rowindex)

end subroutine testcoocsr
!ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

After calling the above subroutine, I got

amtrx=(0, 0, 0, 0, 0, 1)

colindex=(1, 1, 2, 0, 0, 0)

rowindex=(3, 3, 4, 4, 4)

I would very much appreciate your help. Thanks.

Likun

mecej4 · ‎04-01-2015

You have the wrong type for the matrix value arrayss amtrx and aval. They have to be REAL or COMPLEX. Likewise, the name of the routine has to match the type of the matrix -- real(4), real(8), complex(4) or complex(8). Read the MKL documentation before you write the code.

Likun_T_ · ‎04-01-2015

Hi mecej4,

Thanks very much. I've changed the datatype to be consistent with mkl_zcsrcoo. However, I still got problems. Here is my code:

!========================================================

      integer, dimension(:), allocatable :: irn, jcn
      integer, dimension(:), allocatable :: colindex, rowindex
      complex*16, dimension(:), allocatable :: aval, amtrx
      integer job(8), nnz, n, info

      nnz=7
      n=4

   job(1)=2
   job(2)=1
   job(3)=1
   job(5)=nnz
   job(6)=0

   irn=(/1,2,3,3,4,2,4/)
   jcn=(/1,2,1,3,4,4,1/)

   aval(1)=dcmplx(1.d0, 0.d0)
   aval(2)=dcmplx(2.d0, 0.d0)
   aval(3)=dcmplx(3.d0, 0.d0)
   aval(4)=dcmplx(4.d0, 0.d0)
   aval(5)=dcmplx(5.d0, 0.d0)
   aval(6)=dcmplx(6.d0, 0.d0)
   aval(7)=dcmplx(7.d0, 0.d0)

   call mkl_zcsrcoo(job, n, amtrx, colindex, rowindex, nnz, aval, irn, jcn, info)

!==============================================================

The results are not correct. (amtrx has a zero-element, and non-zero imaginary part). Please see below:

amtrx=(dcmplx(1.d0, 0.d0), dcmplx(2.d0, 0.d0), dcmplx(6.d0, 3.d0), dcmplx(0.d0, 0.d0), dcmplx(4.d0, 0.d0), dcmplx(5.d0, 0.d0), dcmplx(7.d0, 0.d0))

colindex=(1,2,4,1,3,1,4)

rowindex=(1,2,4,6,8)

However, if remove the 7-th entry (delete irn(7), jcn(7), aval(7) and set nnz=6), then everything works fine. Your comment is well appreciated.

Thanks,

Likun

mecej4 · ‎04-01-2015

The second code that you posted is incomplete (the allocate statements and the output statements are missing, in particular), so I cannot tell why/if you obtained the incorrect results that you reported.

Likun_T_ · ‎04-01-2015

I am sorry about that. Below is the complete code:

!===================================================

PROGRAM test

      integer, dimension(:), allocatable :: irn, jcn
      integer, dimension(:), allocatable :: colindex, rowindex
      complex*16, dimension(:), allocatable :: aval, amtrx

      integer job(8), nnz, n, info

      nnz=7
      n=4

     allocate(irn(nnz))
     allocate(jcn(nnz))
     allocate(aval(nnz))

     allocate(colindex(nnz))
     allocate(amtrx(nnz))
     allocate(rowindex(n+1))

     !1 0 0 0
     !0 2 0 6
     !3 0 4 0
     !7 0 0 5

   job(1)=2 !from COO to CSR, column indices in CSRare sorted in increasing order within each row.
   job(2)=1 !one-based indexing for the matrix in CSR format is used.
   job(3)=1 !one-based indexing for the matrix in coordinate format is used.
   job(5)=nnz !sets number of the non-zero elements of the matrix A if job(1)=1.
   job(6)=0 !all arrays acsr, ja, ia are filled in for the output storage.

   irn=(/1,2,2,3,3,4,4/)
   jcn=(/1,4,2,3,1,1,4/)

   aval(1)=dcmplx(1.d0, 0.d0)
   aval(2)=dcmplx(6.d0, 0.d0)
   aval(3)=dcmplx(2.d0, 0.d0)
   aval(4)=dcmplx(4.d0, 0.d0)
   aval(5)=dcmplx(3.d0, 0.d0)
   aval(6)=dcmplx(0.d0, -1e-10)
   aval(7)=dcmplx(5.d0, 0.d0)

   call mkl_zcsrcoo(job, n, amtrx, colindex, rowindex, nnz, aval, irn, jcn, info)

   if(info.ne.0) print *, 'conversion from COO to CSR failed'

   deallocate(irn)
     deallocate(jcn)
     deallocate(aval)
     deallocate(colindex)
     deallocate(amtrx)
     deallocate(rowindex)

STOP
END

!=========================================

I tested several other cases, where i found that if irn is sorted in increasing order and jcn is also sorted in increasing order within each row, then the results are correct. However, if jcn is not sorted (as the above code), the results are not correct. I am not sure if there is any restriction for COO form.

Thanks very much for your help.

Likun

mecej4 · ‎04-01-2015

The matrices in #3 and #5 are slightly different. You did not state what results you obtained and what you expected for the matrix in #5. Better focus on one test problem where you think MKL is giving wrong results.

I have used these COO/CSR conversion routines (with real matrices) for a number of years and have never seen an error. COO format data do not need to be in any specific order, but CSR format data do.

Likun_T_ · ‎04-02-2015

Hi mecej4,

Thanks very much for your comment. I tested different datatypes (real, complex and double complex). I found that my code works fine for real and complex, but not for double complex. Below are two codes I testes, using mkl_ccsrcoo and mkl_zcsrcoo, respectively. I kept everything the same in both codes, except the datatype.

Test 1:

!=====================================

PROGRAM test

      integer, dimension(:), allocatable :: irn, jcn
      integer, dimension(:), allocatable :: colindex, rowindex
      complex, dimension(:), allocatable :: aval, amtrx

      integer job(8), nnz, n, info

      nnz=7
      n=4

     allocate(irn(nnz))
     allocate(jcn(nnz))
     allocate(aval(nnz))
     allocate(colindex(nnz))
     allocate(amtrx(nnz))
     allocate(rowindex(n+1))

   job(1)=2
   job(2)=1
   job(3)=1
   job(5)=nnz
   job(6)=0

   irn=(/1,2,2,3,3,4,4/)
   jcn=(/1,4,2,3,1,1,4/)

   aval(1)=cmplx(1,8)
   aval(2)=cmplx(6,0)
   aval(3)=cmplx(2,7)
   aval(4)=cmplx(4,4)
   aval(5)=cmplx(3,2)
   aval(6)=cmplx(0,1)
   aval(7)=cmplx(5,3)

   call mkl_ccsrcoo(job, n, amtrx, colindex, rowindex, nnz, aval, irn, jcn, info)

   if(info.ne.0) print *, 'conversion from COO to CSR failed'

   deallocate(irn)
     deallocate(jcn)
     deallocate(aval)
     deallocate(colindex)
     deallocate(amtrx)
     deallocate(rowindex)

STOP
END

!============================================

The results are

amtrx=(/(1,8), (2,7), (6,0), (3,2), (4,4), (0,1), (5,3)/)

colindex=(/1,2,4,1,3,1,4/)

rowlindex=(/1,2,4,6,8/)

which are in CSR form.

Test 2:

!=====================================

PROGRAM test

      integer, dimension(:), allocatable :: irn, jcn
      integer, dimension(:), allocatable :: colindex, rowindex
      double complex, dimension(:), allocatable :: aval, amtrx

      integer job(8), nnz, n, info

      nnz=7
      n=4

     allocate(irn(nnz))
     allocate(jcn(nnz))
     allocate(aval(nnz))
     allocate(colindex(nnz))
     allocate(amtrx(nnz))
     allocate(rowindex(n+1))

   job(1)=2
   job(2)=1
   job(3)=1
   job(5)=nnz
   job(6)=0

   irn=(/1,2,2,3,3,4,4/)
   jcn=(/1,4,2,3,1,1,4/)

   aval(1)=dcmplx(1,8)
   aval(2)=dcmplx(6,0)
   aval(3)=dcmplx(2,7)
   aval(4)=dcmplx(4,4)
   aval(5)=dcmplx(3,2)
   aval(6)=dcmplx(0,1)
   aval(7)=dcmplx(5,3)

   call mkl_zcsrcoo(job, n, amtrx, colindex, rowindex, nnz, aval, irn, jcn, info)

   if(info.ne.0) print *, 'conversion from COO to CSR failed'

   deallocate(irn)
     deallocate(jcn)
     deallocate(aval)
     deallocate(colindex)
     deallocate(amtrx)
     deallocate(rowindex)

STOP
END

!============================================

The results are

amtrx=(/(1,6), (8,2), (0,7), (4,4), (3,2), (0,1), (5,3)/)

colindex=(/1,2,4,1,3,1,4/)

rowlindex=(/1,2,4,6,8/)

which are not correct.

It seems my definition of 'double complex' is problematic, but I did not see the reason. I would very much appreciate your help. Thanks.

Likun

mecej4 · ‎04-02-2015

You changed the matrix again, without any apparent reason. You have no WRITE statements for the results, so I do not know how you obtained the incorrect results that you reported. Nor have you specified which compiler and MKL versions you used and the compiler options in effect. It is not possible to help you with this problem if you fail to provide clear and complete reports and instructions regarding how to reproduce the problem.

I ran both the codes in #7, after adding WRITE statements to output the CSR representation, and obtained correct results with both the 14.0.4.237 and 15.0.2.179 compilers in 32-bit mode, using the MKL versions that came with those compilers.

Brian_Murphy · ‎09-14-2015

I have found that mkl_zcsrcoo behaves identically as mkl.ccsrcoo. If you feed either one complex*8 input arrays, they produce the same output, which is correct. If you feed complex*16 arrays to the z routine, you get wrong output.