topic Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices in Intel® oneAPI Math Kernel Library

mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Razmo86 — Wed, 18 Aug 2021 19:59:24 GMT

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

Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Gennady_F_Intel — Thu, 19 Aug 2021 04:12:16 GMT

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.

Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Razmo86 — Thu, 19 Aug 2021 12:05:56 GMT

@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

Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Razmo86 — Thu, 19 Aug 2021 18:57:42 GMT

@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

Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Kirill_V_Intel — Thu, 19 Aug 2021 21:47:02 GMT

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

Re: mkl_sparse_sypr for symmetric (Non-Hermitian) complex matrices

Razmo86 — Fri, 20 Aug 2021 13:01:12 GMT

@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

Re:mkl_sparse_sypr doesn't return correct answer for complex matrices

Gennady_F_Intel — Fri, 20 Aug 2021 14:27:47 GMT

Razmo, the official channel to submit the Feature Request to the Intel Online Service Center - https://supporttickets.intel.com/servicecenter?lang=en-US

Re:mkl_sparse_sypr doesn't return correct answer for complex matrices

Gennady_F_Intel — Fri, 01 Oct 2021 03:26:08 GMT

Intel oneAPI Math Kernel Library (oneMKL) version 2021.4 is now available. The MKL Developer Reference has been updated.

Re:mkl_sparse_sypr doesn't return correct answer for complex matrices

Gennady_F_Intel — Fri, 01 Oct 2021 03:26:23 GMT

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.