mkl_sparse_z_syprd, fortran

Bertazzi__Francesco · ‎07-10-2018

It seems that mkl_sparse_z_syprd requires real arguments even though they should be complex. Any advice? Thanks!

Spencer_P_Intel · ‎07-10-2018

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.

Spencer_P_Intel · ‎07-10-2018

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?

Bertazzi__Francesco · ‎07-10-2018

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

Spencer_P_Intel · ‎07-10-2018

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

Bertazzi__Francesco · ‎07-11-2018

It works, thanks!

Cheers, Francesco