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

## MKL-SPBLAS + DSS example : segfault for big matrices

Beginner
1,385 Views

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```

1 Solution
Moderator
1,385 Views

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

4 Replies
Moderator
1,385 Views

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?

Beginner
1,385 Views

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```

Moderator
1,386 Views

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

Beginner
1,385 Views

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