Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

conversion from COO format to CSR

Likun_T_
Beginner
942 Views

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

 

0 Kudos
8 Replies
mecej4
Honored Contributor III
942 Views

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.

0 Kudos
Likun_T_
Beginner
942 Views

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
   

0 Kudos
mecej4
Honored Contributor III
942 Views

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.

0 Kudos
Likun_T_
Beginner
942 Views

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

 

0 Kudos
mecej4
Honored Contributor III
942 Views

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.

0 Kudos
Likun_T_
Beginner
942 Views

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

 

0 Kudos
mecej4
Honored Contributor III
942 Views

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.

0 Kudos
Brian_Murphy
New Contributor II
942 Views

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.

 

0 Kudos
Reply