customized cusparse module

Hello all,

Recently, i wrote my own code to use cusparse module while interfacing with C library.

The code was compiled well but executed in a wrong way.

Can you take a look at my code for what is wrong?

Much appreciate it.

Kevin

    module cusparse_m
      use iso_c_binding

      enum, bind(C) ! cusparseStatus_t
          enumerator :: CUSPARSE_STATUS_SUCCESS=0
          enumerator :: CUSPARSE_STATUS_NOT_INITIALIZED=1
          enumerator :: CUSPARSE_STATUS_ALLOC_FAILED=2
          enumerator :: CUSPARSE_STATUS_INVALID_VALUE=3
          enumerator :: CUSPARSE_STATUS_ARCH_MISMATCH=4
          enumerator :: CUSPARSE_STATUS_MAPPING_ERROR=5
          enumerator :: CUSPARSE_STATUS_EXECUTION_FAILED=6
          enumerator :: CUSPARSE_STATUS_INTERNAL_ERROR=7
          enumerator :: CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED=8
      end enum

      enum, bind(C) ! cusparseMatrixType_t 
         enumerator :: CUSPARSE_MATRIX_TYPE_GENERAL = 0
         enumerator :: CUSPARSE_MATRIX_TYPE_SYMMETRIC = 1
         enumerator :: CUSPARSE_MATRIX_TYPE_HERMITIAN = 2
         enumerator :: CUSPARSE_MATRIX_TYPE_TRIANGULAR = 3
      end enum

      enum, bind(C) ! cusparseIndexBase_t 
         enumerator :: CUSPARSE_INDEX_BASE_ZERO = 0
         enumerator :: CUSPARSE_INDEX_BASE_ONE = 1
      end enum

      enum, bind(C) ! cusparseDirection_t 
         enumerator :: CUSPARSE_DIRECTION_ROW = 0
         enumerator :: CUSPARSE_DIRECTION_COLUMN = 1
      end enum

      enum, bind(C) ! cusparseOperation_t 
         enumerator :: CUSPARSE_OPERATION_NON_TRANSPOSE = 0
         enumerator :: CUSPARSE_OPERATION_TRANSPOSE = 1
         enumerator :: CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE = 2
      end enum

      type cusparseHandle
        type(c_ptr) :: handle
      end type cusparseHandle

      type cusparseMatDescr
        type(c_ptr) :: descr
      end type cusparseMatDescr

      type cusparseSolveAnalysisInfo
        type(c_ptr) :: info
      end type cusparseSolveAnalysisInfo

      interface cusparseCreate
         integer function cusparseCreate(handle) &
                 bind(C,name="cusparseCreate")
            use iso_c_binding
            import cusparseHandle
            implicit none
            type(cusparseHandle) :: handle
         end function cusparseCreate
      end interface

      interface cusparseGetVersion
         integer function cusparseGetVersion(handle, version) &
                 bind(C,name="cusparseGetVersion")
            use iso_c_binding
            import cusparseHandle
            implicit none
            type(cusparseHandle) :: handle
            integer(c_int)       :: version
         end function cusparseGetVersion
      end interface

      interface cusparseCreateMatDescr
         integer function cusparseCreateMatDescr(descr) &
                 bind(C,name="cusparseCreateMatDescr")
            use iso_c_binding
            import cusparseMatDescr
            implicit none
            type(cusparseMatDescr) :: descr
         end function cusparseCreateMatDescr
      end interface

      interface cusparseSetMatType
         integer function cusparseSetMatType(descr, type) &
                 bind(C,name="cusparseSetMatType")
            use iso_c_binding
            import cusparseMatDescr
            implicit none
            type(cusparseMatDescr) :: descr
            integer(c_int) :: type
         end function cusparseSetMatType
      end interface

      interface cusparseSetMatIndexBase
         integer function cusparseSetMatIndexBase(descr, base) &
                 bind(C,name="cusparseSetMatType")
            use iso_c_binding
            import cusparseMatDescr
            implicit none
            type(cusparseMatDescr) :: descr
            integer(c_int) :: base
         end function cusparseSetMatIndexBase
      end interface

      interface cusparseCreateSolveAnalysisInfo
         integer function cusparseCreateSolveAnalysisInfo(info) &
                 bind(C,name="cusparseCreate")
            use iso_c_binding
            import cusparseSolveAnalysisInfo
            implicit none
            type(cusparseSolveAnalysisInfo) :: info
         end function cusparseCreateSolveAnalysisInfo
      end interface

      interface cusparseDestroy
         integer function cusparseDestroy(handle) &
                 bind(C,name="cusparseDestroy")
            use iso_c_binding
            import cusparseHandle
            implicit none
            type(cusparseHandle) :: handle
         end function cusparseDestroy
      end interface

      interface cusparseDnnz
         integer function cusparseDnnz(handle, dirA, m, n, descr, Aval, lda, &
                 nnzPRC, nnzTDHP) bind(C,name="cusparseDnnz")
            use iso_c_binding
            import cusparseHandle
            import cusparseMatDescr
            implicit none
            type(cusparseHandle) :: handle
            type(cusparseMatDescr) :: descr
            integer(c_int) :: dirA, m, n, lda, nnzTDHP
            integer(c_int), device :: nnzPRC(:)
            real(c_double), device :: Aval(:,:)
         end function cusparseDnnz
      end interface

      interface cusparseDdense2csr
         integer function cusparseDdense2csr(handle, m, n, descr, Aval, lda, &
                 nnzPR, csrValA, csrRowPA, csrColIA) bind(C,name="cusparseDdense2csr")
            use iso_c_binding
            import cusparseHandle
            import cusparseMatDescr
            implicit none
            type(cusparseHandle) :: handle
            type(cusparseMatDescr) :: descr
            integer(c_int) :: m, n, lda, nnzTDHP
            integer(c_int), device :: nnzPR(:), csrRowPA(:), csrColIA(:)
            real(c_double), device :: Aval(:,:), csrValA(:)
         end function cusparseDdense2csr
      end interface

      interface cusparseDcsrmm
         integer function cusparseDcsrmm(handle, transA, m, n, k, nnz, alpha, descr, &
                 csrValA, csrRowPA, csrColIA, B, ldb, beta, C, ldc) bind(C,name="cusparseDcsrmm")
            use iso_c_binding
            import cusparseHandle
            import cusparseMatDescr
            implicit none
            type(cusparseHandle) :: handle
            type(cusparseMatDescr) :: descr
            integer(c_int) :: transA, m, n, k, nnz, ldb, ldc
            real(c_double) :: alpha, beta
            integer(c_int), device :: csrRowPA(:), csrColIA(:)
            real(c_double), device :: B(:,:), C(:,:), csrValA(:)
         end function cusparseDcsrmm
      end interface

    end module cusparse_m

! Compile with "pgfortran testLevel3.cuf -lcusparse"    
    program testLevel3
    use cudafor
    use cusparse_m
    
    implicit none

    integer, parameter :: nd=20 ! # rows/cols in dense matrix

    type(cusparseHandle) :: h
    type(cusparseMatDescr) :: descrA
    type(cusparseSolveAnalysisInfo) :: saInfo
    integer :: status, version, mode, i

!   D-data
!   dense
    real(8) :: DAde(nd,nd), DBde(nd,nd), DCde(nd,nd), DmaxErr
    real(8), device :: DAde_d(nd,nd), DBde_d(nd,nd), DCde_d(nd,nd)
!   csr
    real(8) :: csrValDA(nd)
    real(8), device :: csrValDA_d(nd)
    real(8) :: Dalpha, Dbeta
    real(8), device :: Dalpha_d, Dbeta_d

!   integer data common to all data types
    integer :: nnz
    integer :: nnzPerRowA(nd), csrRowPtrA(nd+1), csrColIndA(nd)
    integer, device :: nnzPerRowA_d(nd), csrRowPtrA_d(nd+1), csrColInda_d(nd)

!   Initialization

    status = cusparseCreate(h)
    write(*,*) '1', status
    status = cusparseGetVersion(h, version)
    write(*,*) '2', status
    status = cusparseCreateMatDescr(descrA)
    write(*,*) '3', status
    status = cusparseSetMatType(descrA,CUSPARSE_MATRIX_TYPE_GENERAL)
    write(*,*) '4', status
    status = cusparseSetMatIndexBase(descrA,CUSPARSE_INDEX_BASE_ONE)
    write(*,*) '5', status
    status = cusparseCreateSolveAnalysisInfo(saInfo)
    write(*,*) '6', status

!   Initialize matrix (Identity)
   
    DAde = 0.0
    do i = 1, nd
       DAde(i,i) = 1.0
    end do
    DAde_d = DAde
    call random_number(DBde)
    DBde_d = DBde

!   convert from dense to csr
    status = cusparseDnnz(h, CUSPARSE_DIRECTION_ROW, nd, nd, descrA, &
             DAde_d, nd, nnzPerRowA_d, nnz)
    status = cusparseDdense2csr(h, nd, nd, descrA, DAde_d, nd, nnzPerRowA_d, &
             csrValDA_d, csrRowPtrA_d, csrColIndA_d)

!   csrmm HPM
    Dalpha = 1.0
    Dbeta = 0.0
    status = cusparseDcsrmm(h, CUSPARSE_OPERATION_NON_TRANSPOSE, nd, nd, nd, &
             nnz, Dalpha, descrA, csrValDA_d, csrRowPtrA_d, csrColIndA_d, DBde_d, &
             nd, Dbeta, DCde_d, nd)
    if (status /= CUSPARSE_STATUS_SUCCESS) write (*,*) 'CSRMM Error:', status

    DCde = DCde_d
    DmaxErr = maxval(abs(DCde-DBde))

    status = cusparseDestroy(h)
    write(*,*) 'cusparseDestroy', status, DmaxErr

end program testLevel3

After compiling this code, you can see the following result.

1 0
2 1
3 0
4 3
5 3
6 0
CSRMM Error : 1

Hi Kevin,

I see a few problems where you’re not using the “value” attribute for arguments that need to be passed by value to the cuSparse routines. Also, you need to be passing arrays as assumed-size (A(*)) instead of assumed-shape (A(:,:)). Assumed-shape arrays also pass in the descriptor while assumed-size only pass in the pointer to the data. For multi-dimensional arrays, you also want to use a ignore_tkr directive so that the compiler wont give a syntax error. For example:

      !pgi$ ignore_tkr(r) A

I highly suggest using the cusparse module that we ship with the compilers. cuSparse is a tricky one to get correct.

  • Mat

Hi Mat

It helps a lot.

Thank you.

Merry Christmas.

Kevin

Hi,

I am also trying to use cuSparse with Fortran OpenACC.

In order to avoid all the extra source code, I decided to
write the cuSparse parts in a C file and jsut call the C functions from Fortran with an acc host_data use_device.

I use the iso C bindings to avoid compiler-dependent interface issues.

If you would like a sample of my code, let me know.

I would also be interested in comparing notes on the performance of the cusparse solve calls.