- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I need to compute the symmetric product of three sparse matrices (D=P'AP) using the mkl_sparse_sypr routine, however I am not able to get the correct results. In my case, the P matrix is a general type real matrix, and the matrix A is a symmetric complex matrix. For example, consider the below matrices and their corresponding 3-array CSR format (for matrix A only the upper triangle is stored):
For matrix A, I created the CSR matrix and its descriptor as shown below:
Similarly, for matrix P, it is defined as a general type matrix. Here I declared the REAL matrix P as a complex matrix ( the imaginary part of all the nonzero entries are set to zero):
Finally, the below code is used to perform the symmetric product in a single computation stage. I also export the CSR results in order to validate the entries of the resultant symmetric matrix D.
The resultant matrix D has a correct sparsity pattern, however, in this example, two entries doesn't have a correct value (see highlighted below). The correct values are -4+12i and 2+12i, respectively.
I observed the same issue once I tried to compute the PAP' using the SPARSE_OPERATION_NON_TRANSPOSE operation. The results appear to be correct for the cases where both A and P matrices are real.
Is mkl_sparse_sypr routine only applicable to Hermitian complex matrix A ? Can it be applied to symmetric complex matrix A?
Thanks
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This routine (C:=A*B*opA(A) ) is applicable if A has a general structure, B and C are symmetric or Hermitian matrices. We don't know the issue wrt this routine in the current version of MKL. if you observe some problem - I would recommend give us the reproducer which we could compile and run on our end.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@Gennady_F_Intel Thank you for looking into this issue.
!============================================================================
!
! Consider the matrix A
!
! | 1 -1+3i 0 -3+6i 0 |
! | -1+3i 5 0 0 8 |
! A = | 0 0 4 6 4 |,
! | -3+6i 0 6 7 0 |
! | 0 8 4 0 -5 |
!
! The matrix A is represented in a one-based compressed sparse row (CSR) storage
! scheme with three arrays: Only the upper triangular part is stored
!
! values = ( 1 -1+3i -3+6i 5 8 4 6 4 7 -5 )
! columns = ( 1 2 4 2 5 3 4 5 4 5 )
! rowIndex = ( 1 4 6 9 10 11 )
!
! | 1 0 0 1 0 |
! | 0 0 4 0 0 |
! P = | 0.5 0 0 0 0 |,
! | 0 0 0 1 0 |
! | 0 -6.5 0 0 1 |
!
! values = ( 1 1 4 0.5 1 -6.5 1 )
! columns = ( 1 4 3 1 4 2 5 )
! rowIndex = ( 1 3 4 5 6 8 )
! The test performs the following operations :
! The test performs the following operation :
!
! The code computes D = P'AP using mkl_sparse_sypr. P is a MXM General
sparse matrix
! where A is a symmetric sparse matrix where only the upper triangular part is stored
!
!*******************************************************************************
PROGRAM TestSPBLASProgram
USE MKL_SPBLAS
USE ISO_C_BINDING
IMPLICIT NONE
INTEGER M, N, NNZ,NNZP, i, info,j
! *****************************************************************************
! Sparse representation of the matrix A and P
! *****************************************************************************
INTEGER, ALLOCATABLE :: csrColInd(:), csrRowPtr(:),csrPColInd(:), csrPRowPtr(:)
COMPLEX*16, ALLOCATABLE :: csrVal(:),csrPVal(:)
! Matrix descriptor
TYPE(MATRIX_DESCR) descrA,descrP,descrD ! Sparse matrix descriptor
! CSR matrix representation
TYPE(SPARSE_MATRIX_T) csrA,csrP,csrD ! Structure with sparse matrix
! *****************************************************************************
! Declaration of local variables:
! *****************************************************************************
INTEGER :: nCol,nrowsD, ncolsD
INTEGER(C_INT) :: indexing
TYPE(C_PTR) :: rowsD_start, rowsD_end, colD_indx, Dvalues
INTEGER , POINTER :: rowsD_start_f(:), rowsD_end_f(:), colD_indx_f(:)
COMPLEX*16, POINTER :: Dvalues_f(:)
NNZ = 10
ALLOCATE(csrColInd(NNZ))
ALLOCATE(csrRowPtr(M+1))
ALLOCATE(csrVal(NNZ))
csrVal = (/(1.0D0,0.0D0), (-1.0D0,3.0D0) ,(-3.0D0,6.0D0), (5.0D0,0.0D0) ,(8.0D0,0.0D0) ,(4.0D0,0.0D0) ,(6.0D0,0.0D0), (4.0D0,0.0D0), (7.0D0,0.0D0), (-5.0D0,0.0D0)/)
csrColInd = (/ 1, 2, 4, 2, 5, 3, 4, 5, 4, 5 /)
csrRowPtr = (/1, 4, 6, 9, 10, 11/)
M = 5
N = 5
NNZP = 6
ALLOCATE(csrPColInd(NNZP))
ALLOCATE(csrPRowPtr(M+1))
ALLOCATE(csrPVal(NNZP))
csrPVal = (/(1.0D0,0.0D0),(1.0D0,0.0D0) ,(4.0D0,0.0D0) ,(0.5D0,0.0D0),(1.0D0,0.0D0) ,(-6.5D0,0.0D0) ,(1.0D0,0.0D0)/)
csrPColInd = (/1 ,4 ,3 ,1 ,4 ,2 ,5 /)
csrPRowPtr = (/1 ,3 ,4 ,5 ,6 ,8 /)
print*,'EXAMPLE PROGRAM FOR mkl_sparse_sypr'
print*,'--------------------------------------------------------'
! Create CSR matrix A
i = MKL_SPARSE_Z_CREATE_CSR(csrA,SPARSE_INDEX_BASE_ONE,M,N,csrRowPtr,csrRowPtr(2),csrColInd,csrVal)
! Create matrix descriptor
!descrA % TYPE = SPARSE_MATRIX_TYPE_GENERAL
descrA % TYPE = SPARSE_MATRIX_TYPE_SYMMETRIC
descrA % MODE = SPARSE_FILL_MODE_UPPER
descrA % DIAG = SPARSE_DIAG_NON_UNIT
! Analyze sparse matrix; chose proper kernels and workload balancing strategy
info = MKL_SPARSE_OPTIMIZE(csrA)
! Create CSR matrix P
i = MKL_SPARSE_Z_CREATE_CSR(csrP,SPARSE_INDEX_BASE_ONE,M,N,csrPRowPtr,csrPRowPtr(2),csrPColInd,csrPVal)
! Create matrix descriptor
descrP % TYPE = SPARSE_MATRIX_TYPE_GENERAL
! Analyze sparse matrix; chose proper kernels and workload balancing strategy
info = MKL_SPARSE_OPTIMIZE(csrP)
! Single Stage Computations
info = mkl_sparse_sypr (SPARSE_OPERATION_CONJUGATE_TRANSPOSE, csrP, csrA, descrA, csrD, SPARSE_STAGE_FULL_MULT)
! export the matrix D in CSR format matrix from the internal representation
info = mkl_sparse_z_export_csr (csrD, indexing, nrowsD, ncolsD, rowsD_start, rowsD_end, colD_indx, Dvalues)
info = mkl_sparse_order(csrD)
! Converting C into Fortran pointers
call C_F_POINTER(rowsD_start, rowsD_start_f, [nrowsD])
call C_F_POINTER(rowsD_end , rowsD_end_f , [nrowsD])
call C_F_POINTER(colD_indx , colD_indx_f , [rowsD_end_f(nrowsD)-indexing])
call C_F_POINTER(Dvalues , Dvalues_f , [rowsD_end_f(nrowsD)-indexing])
! Release internal representation of CSR matrix
info = MKL_SPARSE_DESTROY(csrA)
info = MKL_SPARSE_DESTROY(csrP)
info = MKL_SPARSE_DESTROY(csrD)
END PROGRAM TestSPBLASProgram
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@Gennady_F_Intel Thanks for looking into this issue. Here is the test program:
!===============================================================================
!
! Consider the matrix A
!
! | 1 -1+3i 0 -3+6i 0 |
! | -1+3i 5 0 0 8 |
! A = | 0 0 4 6 4 |,
! | -3+6i 0 6 7 0 |
! | 0 8 4 0 -5 |
!
! The matrix A is represented in a one-based compressed sparse row (CSR) storage
! scheme with three arrays: Only the upper triangular part is stored
!
! values = ( 1 -1+3i -3+6i 5 8 4 6 4 7 -5 )
! columns = ( 1 2 4 2 5 3 4 5 4 5 )
! rowIndex = ( 1 4 6 9 10 11 )
!
! | 1 0 0 1 0 |
! | 0 0 4 0 0 |
! P = | 0.5 0 0 0 0 |,
! | 0 0 0 1 0 |
! | 0 -6.5 0 0 1 |
!
! values = ( 1 1 4 0.5 1 -6.5 1 )
! columns = ( 1 4 3 1 4 2 5 )
! rowIndex = ( 1 3 4 5 6 8 )
! The test performs the following operations :
!
! The test performs the following operations :
!
! The code computes D = P'AP using mkl_sparse_sypr. P is a MXM General sparse matrix
! where A is a symmetric sparse matrix where only the upper triangular part is stored
!
!*******************************************************************************
PROGRAM TestSPBLASProgram
USE MKL_SPBLAS
USE ISO_C_BINDING
IMPLICIT NONE
INTEGER M, N, NNZ,NNZP, i, info,j
! *****************************************************************************
! Sparse representation of the matrix A and P
! *****************************************************************************
INTEGER, ALLOCATABLE :: csrColInd(:), csrRowPtr(:),csrPColInd(:), csrPRowPtr(:)
COMPLEX*16, ALLOCATABLE :: csrVal(:),csrPVal(:)
! Matrix descriptor
TYPE(MATRIX_DESCR) descrA,descrP,descrD ! Sparse matrix descriptor
! CSR matrix representation
TYPE(SPARSE_MATRIX_T) csrA,csrP,csrD ! Structure with sparse matrix
! *****************************************************************************
! Declaration of local variables:
! *****************************************************************************
INTEGER :: nCol,nrowsD, ncolsD
INTEGER(C_INT) :: indexing
TYPE(C_PTR) :: rowsD_start, rowsD_end, colD_indx, Dvalues
INTEGER , POINTER :: rowsD_start_f(:), rowsD_end_f(:), colD_indx_f(:)
COMPLEX*16, POINTER :: Dvalues_f(:)
NNZ = 10
ALLOCATE(csrColInd(NNZ))
ALLOCATE(csrRowPtr(M+1))
ALLOCATE(csrVal(NNZ))
csrVal = (/(1.0D0,0.0D0), (-1.0D0,3.0D0) ,(-3.0D0,6.0D0), (5.0D0,0.0D0) ,(8.0D0,0.0D0) ,(4.0D0,0.0D0) ,(6.0D0,0.0D0), (4.0D0,0.0D0), (7.0D0,0.0D0), (-5.0D0,0.0D0)/)
csrColInd = (/ 1, 2, 4, 2, 5, 3, 4, 5, 4, 5 /)
csrRowPtr = (/1, 4, 6, 9, 10, 11/)
M = 5
N = 5
NNZP = 6
ALLOCATE(csrPColInd(NNZP))
ALLOCATE(csrPRowPtr(M+1))
ALLOCATE(csrPVal(NNZP))
csrPVal = (/(1.0D0,0.0D0),(1.0D0,0.0D0) ,(4.0D0,0.0D0) ,(0.5D0,0.0D0),(1.0D0,0.0D0) ,(-6.5D0,0.0D0) ,(1.0D0,0.0D0)/)
csrPColInd = (/1 ,4 ,3 ,1 ,4 ,2 ,5 /)
csrPRowPtr = (/1 ,3 ,4 ,5 ,6 ,8 /)
print*,'EXAMPLE PROGRAM FOR mkl_sparse_sypr'
print*,'--------------------------------------------------------'
! Create CSR matrix A
i = MKL_SPARSE_Z_CREATE_CSR(csrA,SPARSE_INDEX_BASE_ONE,M,N,csrRowPtr,csrRowPtr(2),csrColInd,csrVal)
! Create matrix descriptor
!descrA % TYPE = SPARSE_MATRIX_TYPE_GENERAL
descrA % TYPE = SPARSE_MATRIX_TYPE_SYMMETRIC
descrA % MODE = SPARSE_FILL_MODE_UPPER
descrA % DIAG = SPARSE_DIAG_NON_UNIT
! Analyze sparse matrix; chose proper kernels and workload balancing strategy
info = MKL_SPARSE_OPTIMIZE(csrA)
! Create CSR matrix P
i = MKL_SPARSE_Z_CREATE_CSR(csrP,SPARSE_INDEX_BASE_ONE,M,N,csrPRowPtr,csrPRowPtr(2),csrPColInd,csrPVal)
! Create matrix descriptor
descrP % TYPE = SPARSE_MATRIX_TYPE_GENERAL
! Analyze sparse matrix; chose proper kernels and workload balancing strategy
info = MKL_SPARSE_OPTIMIZE(csrP)
! Single Stage Computations
info = mkl_sparse_sypr (SPARSE_OPERATION_CONJUGATE_TRANSPOSE, csrP, csrA, descrA, csrD, SPARSE_STAGE_FULL_MULT)
! export the matrix D in CSR format matrix from the internal representation
info = mkl_sparse_z_export_csr (csrD, indexing, nrowsD, ncolsD, rowsD_start, rowsD_end, colD_indx, Dvalues)
info = mkl_sparse_order(csrD)
! Converting C into Fortran pointers
call C_F_POINTER(rowsD_start, rowsD_start_f, [nrowsD])
call C_F_POINTER(rowsD_end , rowsD_end_f , [nrowsD])
call C_F_POINTER(colD_indx , colD_indx_f , [rowsD_end_f(nrowsD)-indexing])
call C_F_POINTER(Dvalues , Dvalues_f , [rowsD_end_f(nrowsD)-indexing])
! Release internal representation of CSR matrix
info = MKL_SPARSE_DESTROY(csrA)
info = MKL_SPARSE_DESTROY(csrP)
info = MKL_SPARSE_DESTROY(csrD)
END PROGRAM TestSPBLASProgram
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi!
I'm pretty sure the reason for what you see is a misleading wording in documentation. Routine mkl_sparse_sypr supports symmetric matrices A for real data types and hermitian (only) A for complex case [while the docs say "B and C are symmetric or Hermitian matrices" as if symmetric works for complex matrices].
E.g., a quick implicit proof: if you change
descrA % TYPE = SPARSE_MATRIX_TYPE_SYMMETRIC
to
descrA % TYPE = SPARSE_MATRIX_TYPE_HERMITIAN
you'll get same result. I haven't checked but I believe that is the correct result if matrix A is treated as Hermitian [let me know if you think this is incorrect].
If you need this functionality, our TCE can help you submit a FR, so that we can extend mkl-sparse_sypr to work for your case (should be straightforward).
Best,
Kirill
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
@Kirill_V_Intel Thanks for your answer. Correct, I have already verified that the routine is working correctly for real symmetric and Hermitian complex matrices. It would be great if you could extend mkl-sparse_sypr to symmetric complex matrices as well. Would you please let me know how I should process with the FR submission?
Regards
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Razmo, the official channel to submit the Feature Request to the Intel Online Service Center - https://supporttickets.intel.com/servicecenter?lang=en-US
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Intel oneAPI Math Kernel Library (oneMKL) version 2021.4 is now available. The MKL Developer Reference has been updated.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This issue has been resolved and we will no longer respond to this thread. If you require additional assistance from Intel, please start a new thread. Any further interaction in this thread will be considered community only.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page