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

## mkl_sparse_z_syprd, fortran

Beginner
504 Views

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

5 Replies
Employee
504 Views

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.

Employee
504 Views

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?

Beginner
504 Views

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

Employee
504 Views

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:

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

Beginner
504 Views

It works, thanks!

Cheers, Francesco