Solved: Thank you, now it almost

Tibaldi__Alberto · ‎01-24-2020

Hi,

I am a newbie with MKL.

I am trying to combine SPBLAS and DSS on a matrix (which is basically the mass matrix of a 1D FEM solver) which I can handle better in COO format. Therefore, what I've done is: store it with COO format in SPBLAS, convert such internal representation to CSR, perform a matrix-vector product in the internal representation, convert to the three-array variation of CSR (three vectors). With this, I can generate a DSS handle, process it, and finally solve a system.

Before starting with the hard stuff, I decided to solve a canonical system such that its solution is equal to all ones, so I can check the solution correctness (to verify that all steps are performed correctly).

It appears that everything works, but just for "small" matrices. For instance, nn=100 works, nn=10 works, nn=5000 doesn't, and segfaults.

I have therefore some questions.

1) Is what I'm trying to do the best way of solving a sparse linear system starting from a COO matrix?

2) Any clue about why it doesn't it work with bigger (but not so big.....) matrices?

3) I noticed that, recently, a QR sparse solver working directly with the internal SPBLAS representation has been developed. I will try it ASAP, but I would have preferred using DSS/PARDISO ... And I don't know which is the best way, see question 1.
Is there any chance to make the SPBLAS internal representation directly "digestible" by DSS or PARDISO? Now, or, in the future?

Please find attached the code I've cooked from the MKL examples...

Thank you!

Alberto

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This code assembles the typical mass matrix of a 1-D FEM simulator in COO - SPBLAS, converts
! it to CSR, computes the known term such that Mmat*xref = 1, with MKL_SPARSE routines,
! exports the matrix in the 3-array variation of CSR, assembles it in DSS, and solves the
! system achieving all the solution vector components equal to 1.
!
! Alberto Tibaldi, 2020-01-23
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Modules to be included
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
INCLUDE 'mkl_spblas.f90' ! include sparse BLAS module
INCLUDE 'mkl_dss.f90'    ! include DSS module
!
PROGRAM TEST_MKL_SPBLAS_DSS
!
    ! Used modules
    USE ISO_C_BINDING
    USE MKL_SPBLAS
    USE MKL_DSS
!
    IMPLICIT NONE
!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Variable declarations
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Parameters
    REAL(KIND=8), PARAMETER                  :: pi = 3.141592653589793D0      ! pi 
    INTEGER(KIND=4), PARAMETER               :: nn = 100                        ! number of mesh nodes
!
    ! Main program variables
    INTEGER(KIND=4)                          :: ne = nn - 1 ! number of rows/columns of the matrix
    INTEGER(KIND=4)                          :: nnz, perm(1), nrhs
    INTEGER(KIND=4)                          :: ii, rows, cols
!
    REAL(KIND=8)                             :: zvet(nn), Le(nn-1) ! mesh vectors
    !
    ! COO matrix variables
    INTEGER(KIND=4)                          :: in1(nn-1), in2(nn-1), inr(4*(nn-1)), inc(4*(nn-1)) ! indexes for COO representation
    REAL(KIND=8)                             :: Mcoovalues(4*(nn-1)) ! 
!
    ! CSR matrix variables
    REAL(KIND=8)                             :: xref(nn), rhs(nn), xsol(nn)
    INTEGER(KIND=4), ALLOCATABLE             :: rowIndex(:), rows_end(:), columns(:)
    REAL(KIND=8), ALLOCATABLE                :: values(:)
    
    REAL(KIND=8)                             :: step = pi
!
    ! MKL-related variables
    INTEGER(KIND=4)                          :: mkl_indexing, mkl_out
    REAL(KIND=8)                             :: alpha=1.0D0, beta=0.0D0 ! for matrix-vector product
    TYPE(C_PTR)                              :: rowIndex_C, rows_end_C, columns_C, values_C ! C pointers
    INTEGER, POINTER, DIMENSION(:)           :: rowIndex_F, rows_end_F, columns_F ! FORTRAN pointers to integers
    REAL(KIND=8), POINTER, DIMENSION(:)      :: values_f ! FORTRAN pointer to real
!
    TYPE(SPARSE_MATRIX_T)                    :: MCOO, MCSR ! sparse matrix handles for SPBLAS
    TYPE(MATRIX_DESCR)                       :: descrM     ! sparse matrix descriptor
!
    TYPE(MKL_DSS_HANDLE)                     :: syshandle  ! sparse matrix handle for DSS
!
    descrM%TYPE = SPARSE_MATRIX_TYPE_GENERAL
!
    zvet = step/(nn-1)*(/ (ii, ii=0,nn-1) /)
    in1 = (/ (ii, ii=1,nn-1) /)
    in2 = (/ (ii, ii=2,nn)   /)
    Le = zvet(in2) - zvet(in1)
!
!   Create COO 
    Mcoovalues = (/ 2*Le, 1*Le, 1*Le, 2*Le /)/6
    inr        = (/  in1,  in1,  in2,  in2 /)
    inc        = (/  in1,  in2,  in1,  in2 /)
!
    ! Create COO representation: the three vectors
    nnz = 4*(nn-1) ! overestimation - it will be re-computed later by MKL routines
    mkl_out = MKL_SPARSE_D_CREATE_COO(MCOO,SPARSE_INDEX_BASE_ONE,nn,nn,nnz,inr,inc,Mcoovalues)
    ! Convert MKL-COO into MKL-CSR
    mkl_out = MKL_SPARSE_CONVERT_CSR(MCOO, SPARSE_OPERATION_NON_TRANSPOSE, MCSR)
    ! Optimize MKL-CSR matrix
    mkl_out = MKL_SPARSE_OPTIMIZE(MCSR)
!
    ! Compute matrix-vector product
    xref = 1.0D0
    rhs = 0.0D0
    mkl_out = MKL_SPARSE_D_MV(SPARSE_OPERATION_NON_TRANSPOSE,alpha,MCSR,descrM,xref,beta,rhs)
    ! WRITE(*,*) rhs
    !
    ! Export to CSR: 4-array version (rows_start = rowIndex, rows_end, col_indx, values)
    !                outputs are C pointers
    mkl_out = MKL_SPARSE_D_EXPORT_CSR(MCSR,mkl_indexing,rows,cols,rowIndex_C,rows_end_C,columns_C,values_C)
!
    ! Deallocate MKL matrices
    mkl_out = MKL_SPARSE_DESTROY(MCSR)
    mkl_out = MKL_SPARSE_DESTROY(MCOO)
!
    ! Number of actually stored values in 3-array variation CSR
!
    ! Cast C pointers that came from export routine to Fortran pointers
    CALL C_F_POINTER(rowIndex_C, rowIndex_F, [rows+1])
    CALL C_F_POINTER(rows_end_C, rows_end_F, [rows])
    CALL C_F_POINTER(columns_C, columns_F, [nnz])
    CALL C_F_POINTER(values_C, values_F, [nnz])
!
    nnz = rowIndex_F(rows+1)-mkl_indexing       
    WRITE(*,*) 'Effective number of elements: ', nnz
!
    ALLOCATE(rowIndex(rows+1),rows_end(rows),columns(nnz),values(nnz))
    rowIndex = rowIndex_F
    columns = columns_F
    values = values_F
!
    ! WRITE(*,*) 'rowIndex: ', rowIndex
    ! WRITE(*,*) 'columns: ', columns
    ! WRITE(*,*) 'values: ', values
!
    ! DSS: input parameters
    perm(1) = 0
    nrhs = 1
!
    mkl_out = DSS_CREATE(syshandle, MKL_DSS_DEFAULTS)
    mkl_out = DSS_DEFINE_STRUCTURE(syshandle, MKL_DSS_NON_SYMMETRIC, rowIndex, nn, nn, columns, nnz)
!
    mkl_out = DSS_REORDER(syshandle, MKL_DSS_DEFAULTS, perm)
    mkl_out = DSS_FACTOR_REAL(syshandle, MKL_DSS_DEFAULTS, values)
!
    mkl_out = DSS_SOLVE_REAL(syshandle, MKL_DSS_DEFAULTS, rhs, nRhs, xsol)
!
    ! Print solution on screen
    ! WRITE(*,*) xsol
!
    ! Delete DSS handle
    mkl_out = DSS_DELETE(syshandle, MKL_DSS_DEFAULTS)
!    
    ! Deallocate ALLOCATABLE variables
    DEALLOCATE(rowIndex, rows_end, columns, values)
!
    WRITE(*,*) 'Norm 2 of the solution error: ', NORM2(xsol-1)    

END PROGRAM TEST_MKL_SPBLAS_DSS

Gennady_F_Intel · ‎01-26-2020

please try to destroy the matrixes handles at the very last stage... and check the problem

View solution in original post

Gennady_F_Intel · ‎01-24-2020

it looks like a bug in DSS API. We will check. Which version of mkl do you use? Did you try to use Pardiso API instead of DSS?

Tibaldi__Alberto · ‎01-24-2020

Gennady,

Thank you very much for your answer.

I'm using "Intel(R) Math Kernel Library Version 2020.0.0 Product Build 20191122 for Intel(R) 64 architecture applications".

I've run another test, and now I don't believe that the problem is in DSS, but rather in MKL_SPARSE_D_EXPORT_CSR. In fact, it segfaults even without all the DSS part.

If you run the code at the end of this post, you can see that there are some strange indexes in the columns vector (the last ones are: 2147483647 4998 4999 5000 2147483647 4999 5000).

A piece of information that I've not provided in my original post is: I'm compiling this code with the string

gfortran prova.f90 -fopenmp -I${MKLROOT}/include -L${MKLROOT}/lib/intel64 -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -lm -ldl -fbounds-check -m64 -O2 -o prova

Please find the code below. If you uncomment line 117, it segfaults.

Thank you again for your attention,

Alberto

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This code assembles the typical mass matrix of a 1-D FEM simulator in COO - SPBLAS, converts
! it to CSR, computes the known term such that Mmat*xref = 1, with MKL_SPARSE routines,
! exports the matrix in the 3-array variation of CSR, assembles it in DSS, and solves the
! system achieving all the solution vector components equal to 1.
!
! Alberto Tibaldi, 2020-01-23
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Modules to be included
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
INCLUDE 'mkl_spblas.f90' ! include sparse BLAS module
INCLUDE 'mkl_dss.f90'    ! include DSS module
!
PROGRAM TEST_MKL_SPBLAS_DSS
!
    ! Used modules
    USE ISO_C_BINDING
    USE MKL_SPBLAS
    USE MKL_DSS
!
    IMPLICIT NONE
!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Variable declarations
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Parameters
    REAL(KIND=8), PARAMETER                  :: pi = 3.141592653589793D0      ! pi 
    INTEGER(KIND=4), PARAMETER               :: nn = 5000                        ! number of mesh nodes
!
    ! Main program variables
    INTEGER(KIND=4)                          :: ne = nn - 1 ! number of rows/columns of the matrix
    INTEGER(KIND=4)                          :: nnz, perm(1), nrhs
    INTEGER(KIND=4)                          :: ii, rows, cols
!
    REAL(KIND=8)                             :: zvet(nn), Le(nn-1) ! mesh vectors
    !
    ! COO matrix variables
    INTEGER(KIND=4)                          :: in1(nn-1), in2(nn-1), inr(4*(nn-1)), inc(4*(nn-1)) ! indexes for COO representation
    REAL(KIND=8)                             :: Mcoovalues(4*(nn-1)) ! 
!
    ! CSR matrix variables
    REAL(KIND=8)                             :: xref(nn), rhs(nn), xsol(nn)
    INTEGER(KIND=4), ALLOCATABLE             :: rowIndex(:), rows_end(:), columns(:)
    REAL(KIND=8), ALLOCATABLE                :: values(:)
    
    REAL(KIND=8)                             :: step = pi
!
    ! MKL-related variables
    CHARACTER(LEN=200)                       :: mkl_str  ! useful to print MKL version
    INTEGER(KIND=4)                          :: mkl_indexing, mkl_out
    REAL(KIND=8)                             :: alpha=1.0D0, beta=0.0D0 ! for matrix-vector product
    TYPE(C_PTR)                              :: rowIndex_C, rows_end_C, columns_C, values_C ! C pointers
    INTEGER, POINTER, DIMENSION(:)           :: rowIndex_F, rows_end_F, columns_F ! FORTRAN pointers to integers
    REAL(KIND=8), POINTER, DIMENSION(:)      :: values_f ! FORTRAN pointer to real
!
    TYPE(SPARSE_MATRIX_T)                    :: MCOO, MCSR ! sparse matrix handles for SPBLAS
    TYPE(MATRIX_DESCR)                       :: descrM     ! sparse matrix descriptor
!
    TYPE(MKL_DSS_HANDLE)                     :: syshandle  ! sparse matrix handle for DSS
!
    descrM%TYPE = SPARSE_MATRIX_TYPE_GENERAL
!
    WRITE(*,*) "*************************************************************************************************************** "
    CALL MKL_GET_VERSION_STRING (mkl_str)
    WRITE(*,'(a)') mkl_str
    WRITE(*,*) "*************************************************************************************************************** "
    WRITE(*,*) 
!    
    zvet = step/(nn-1)*(/ (ii, ii=0,nn-1) /)
    in1 = (/ (ii, ii=1,nn-1) /)
    in2 = (/ (ii, ii=2,nn)   /)
    Le = zvet(in2) - zvet(in1)
!
!   Create COO 
    Mcoovalues = (/ 2*Le, 1*Le, 1*Le, 2*Le /)/6
    inr        = (/  in1,  in1,  in2,  in2 /)
    inc        = (/  in1,  in2,  in1,  in2 /)
!
    ! Create COO representation: the three vectors
    nnz = 4*(nn-1) ! overestimation - it will be re-computed later by MKL routines
    mkl_out = MKL_SPARSE_D_CREATE_COO(MCOO,SPARSE_INDEX_BASE_ONE,nn,nn,nnz,inr,inc,Mcoovalues)
    ! Convert MKL-COO into MKL-CSR
    mkl_out = MKL_SPARSE_CONVERT_CSR(MCOO, SPARSE_OPERATION_NON_TRANSPOSE, MCSR)
    ! Optimize MKL-CSR matrix
    mkl_out = MKL_SPARSE_OPTIMIZE(MCSR)
!
    ! Compute matrix-vector product
    xref = 1.0D0
    rhs = 0.0D0
    mkl_out = MKL_SPARSE_D_MV(SPARSE_OPERATION_NON_TRANSPOSE,alpha,MCSR,descrM,xref,beta,rhs)
    ! WRITE(*,*) rhs
    !
    ! Export to CSR: 4-array version (rows_start = rowIndex, rows_end, col_indx, values)
    !                outputs are C pointers
    mkl_out = MKL_SPARSE_D_EXPORT_CSR(MCSR,mkl_indexing,rows,cols,rowIndex_C,rows_end_C,columns_C,values_C)
!
    ! Deallocate MKL matrices
    mkl_out = MKL_SPARSE_DESTROY(MCSR)
    mkl_out = MKL_SPARSE_DESTROY(MCOO)
!
    ! Number of actually stored values in 3-array variation CSR
!
    ! Cast C pointers that came from export routine to Fortran pointers
    CALL C_F_POINTER(rowIndex_C, rowIndex_F, [rows+1])
    CALL C_F_POINTER(rows_end_C, rows_end_F, [rows])
    CALL C_F_POINTER(columns_C, columns_F, [nnz])
    !CALL C_F_POINTER(values_C, values_F, [nnz])
!
    nnz = rowIndex_F(rows+1)-mkl_indexing       
    WRITE(*,*) 'Effective number of elements: ', nnz
!
    ALLOCATE(rowIndex(rows+1),rows_end(rows),columns(nnz),values(nnz))
    rowIndex = rowIndex_F
    columns = columns_F
    ! values = values_F
!
    ! WRITE(*,*) 'rowIndex: ', rowIndex
    WRITE(*,*) 'columns: ', columns
    ! WRITE(*,*) 'values: ', values
!
END PROGRAM TEST_MKL_SPBLAS_DSS

Gennady_F_Intel · ‎01-26-2020

please try to destroy the matrixes handles at the very last stage... and check the problem

Tibaldi__Alberto · ‎01-26-2020

Thank you, now it seems to work. Did I mentioned that I'm a newbie? :-)

The code still segfaults for bigger nn values, but it appears to be a problem of FORTRAN and not MKL. I will investigate it...

Thank you very much for your support!

Alberto

MKL-SPBLAS + DSS example : segfault for big matrices