- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
It seems that mkl_sparse_z_syprd requires real arguments even though they should be complex. Any advice? Thanks!
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Now, there are some caveats that I realize might be a little more confusing so I will see if I can clear them up:
SYPR does the following operations (with opA = non transpose):
C = beta*C + alpha*A*B*(A^H)
Now this operation has been designed to take in a Hermitian matrix B and C and return a Hermitian matrix C. In particular a necessary (not sufficient) condition to be Hermitian is that it must be real valued on diagonal since B_ii = (B^H)_ii = conj(B_ii). (no imaginary parts)
Now given a hermitian matrix B, if we want the resulting C matrix to also be Hermitian we must only use real valued alpha and beta to preserve this same real diagonal property even though technically we can use any complex alpha and beta values.
As a side note, since we do expect a hermitian matrix B and C, we do not use the full general data. We only use the upper triangular part of B in the computation and return only the upper triangular part of C (the lower part can easily be referenced as conjugate transpose of the upper). This allows us to avoid redundant computations and speed up the operation.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I'm not sure I understand the issue.
The header (in C from mkl_spblas.h) is
sparse_status_t mkl_sparse_z_syprd ( const sparse_operation_t op, const sparse_matrix_t A, const MKL_Complex16 *B, const sparse_layout_t layoutB, const MKL_INT ldb, const MKL_Complex16 alpha, const MKL_Complex16 beta, MKL_Complex16 *C, const sparse_layout_t layoutC, const MKL_INT ldc );
and in fortran from mkl_spblas.f90 is
FUNCTION MKL_SPARSE_Z_SYPRD(op,A,B,layoutB,ldb,alpha,beta,C,layoutC,ldc) & ! BIND(C, name='MKL_SPARSE_Z_SYPRD') USE, INTRINSIC :: ISO_C_BINDING , ONLY : C_INT, C_DOUBLE_COMPLEX IMPORT SPARSE_MATRIX_T INTEGER(C_INT) , INTENT(IN) :: op TYPE(SPARSE_MATRIX_T) , INTENT(IN) :: A REAL(C_DOUBLE_COMPLEX) , INTENT(IN) , DIMENSION(*) :: B INTEGER(C_INT) , INTENT(IN) :: layoutB INTEGER(C_INT) , INTENT(IN) :: ldb REAL(C_DOUBLE_COMPLEX) , INTENT(IN) :: alpha REAL(C_DOUBLE_COMPLEX) , INTENT(IN) :: beta REAL(C_DOUBLE_COMPLEX) , INTENT(OUT), DIMENSION(*) :: C INTEGER(C_INT) , INTENT(IN) :: layoutC INTEGER(C_INT) , INTENT(IN) :: ldc INTEGER(C_INT) MKL_SPARSE_Z_SYPRD END FUNCTION
which definitely uses complex valued arguments for alpha, beta, B input and C output. A should also be created as a complex sparse matrix using mkl_sparse_z_create_csr()... or bsr if that is what you are using. Can you clarify what issue you are observing?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thanks for your help!
I think the problem is that I have defined the complex matrix B as
COMPLEX(KIND=8), ALLOCATABLE :: B(:,:)
while in mkl_spblas.f90, B appears as
REAL(C_DOUBLE_COMPLEX) , INTENT(IN) , DIMENSION(*) :: B
It's not clear to me how I should define a dense complex matrix. Could you comment on this?
Thanks, Francesco
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Francesco,
I see that I was a little blind earlier. You are absolutely right. Your matrix is not compatible with this interface. What you can do is the following:
Add this new header definition to the top of your local code:
FUNCTION MKL_SPARSE_Z_SYPRD_F(op,A,B,layoutB,ldb,alpha,beta,C,layoutC,ldc) & ! BIND(C, name='MKL_SPARSE_Z_SYPRD') USE, INTRINSIC :: ISO_C_BINDING , ONLY : C_INT, C_DOUBLE_COMPLEX IMPORT SPARSE_MATRIX_T INTEGER(C_INT) , INTENT(IN) :: op TYPE(SPARSE_MATRIX_T) , INTENT(IN) :: A COMPLEX(C_DOUBLE_COMPLEX) , INTENT(IN) , DIMENSION(*) :: B INTEGER(C_INT) , INTENT(IN) :: layoutB INTEGER(C_INT) , INTENT(IN) :: ldb COMPLEX(C_DOUBLE_COMPLEX) , INTENT(IN) :: alpha COMPLEX(C_DOUBLE_COMPLEX) , INTENT(IN) :: beta COMPLEX(C_DOUBLE_COMPLEX) , INTENT(OUT), DIMENSION(*) :: C INTEGER(C_INT) , INTENT(IN) :: layoutC INTEGER(C_INT) , INTENT(IN) :: ldc INTEGER(C_INT) MKL_SPARSE_Z_SYPRD_F END FUNCTION
And then use MKL_SPARSE_Z_SYPRD_F() instead of MKL_SPARSE_Z_SYPRD() in your code... This should work for your matrix definition.
Best,
Spencer
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
It works, thanks!
Cheers, Francesco

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