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

MKL_SPARSE_D_TRSM ............ no result when called from fortran

may_ka
Beginner
461 Views

Hi

the program below implements a "csr" class which holds a sparse csr matrix copied from the mkl manual. This matrix is used for Inspector-Executor routines MKL_SPARSE_D_CREATE_CSR, MKL_SPARSE_COPY and MKL_SPARSE_D_TRSM. All routines finish without error. However, MKL_SPARSE_D_TRSM does not seem to solve for anything. Am I doing something wrong here??

include "mkl_spblas.f90"
Module Data_Kind
  Implicit None
  Integer, Parameter :: Large=Selected_Int_Kind(12)
  Integer, Parameter :: Medium=Selected_Int_Kind(8)
  Integer(Medium), Parameter :: Double=Selected_Real_Kind(15,100)
End Module Data_Kind
Module Mod_csr
  use data_kind
  implicit None
  Type :: csr
    integer(large), allocatable :: nrows, ncols
    integer(large), allocatable :: rowpos(:), colpos(:),&
      & pointerB(:), pointerE(:)
    real(double), allocatable :: values(:)
  contains
    Procedure :: iii => Subiii
  end type csr
contains
  Subroutine Subiii(this)
    Implicit None
    Class(csr), intent(inout) :: this
    !!2017 mkl manual page 3233
    allocate(this%nrows,this%ncols,source=5)
    allocate(this%rowpos(6),source=(/1,4,6,9,12,14/))
    allocate(this%colpos(13),source=(/1,2,4,1,2,3,4,5,1,3,4,2,5/))
    allocate(this%pointerB(5),source=(/1,4,6,9,12/))
    allocate(this%pointerE(5),source=(/4,6,9,12,14/))
    allocate(this%values(13),&
      &source=(/1.0D0,-1.0D0,-3.0D0,-2.0D0,5.0D0,&
      &4.0D0,6.0D0,4.0D0,-4.0D0,2.0D0,7.0D0,8.0D0,-5.0D0/))
  end Subroutine Subiii
End Module Mod_Csr
Program Test
  use data_kind
  USE IFPORT
  use Mod_csr, only: csr
  use MKL_SPBLAS
  USE, INTRINSIC :: ISO_C_BINDING
  Implicit none
  Type(csr) :: tscsr
  integer(c_int) :: isstat=0
  Type(Sparse_Matrix_T) :: handle, handle1
  Type(MATRIX_DESCR) :: descr
  real(double), allocatable :: in(:,:), out(:,:)
  outer:block
    descr%TYPE=SPARSE_MATRIX_TYPE_GENERAL
    descr%DIAG=SPARSE_DIAG_NON_UNIT
    call tscsr%iii()
    isstat=MKL_SPARSE_D_CREATE_CSR(&
      &handle,&
      &SPARSE_INDEX_BASE_ONE,&
      &tscsr%nrows,&
      &tscsr%ncols,&
      &tscsr%pointerB,&
      &tscsr%pointerE,&
      &tscsr%colpos,&
      &tscsr%values&
      &)
    if(isstat/=0) Then
      write(*,*) "error 1 ",isstat;exit outer
    End if
    isstat=mkl_sparse_copy(handle,descr,handle1)
    if(isstat/=0) Then
      write(*,*) "error 2 ",isstat;exit outer
    End if
    allocate(in(tscsr%nrows,2),out(tscsr%nrows,2))
    in=1.0;out=0.0
    isstat=MKL_SPARSE_D_TRSM(&
      &SPARSE_OPERATION_NON_TRANSPOSE,&
      &1.0_double,&
      &handle1,&
      &descr,&
      &SPARSE_LAYOUT_COLUMN_MAJOR,&
      &in,&
      &2,&
      &size(in,1),&
      &out,&
      &size(out,1)&
      &)
    if(isstat/=0) Then
      write(*,*) "error 3 ",isstat;exit outer
    end if
    write(*,*) maxval(in), maxval(out), minval(out)
  end block outer
End Program Test

Matrix "out" is supposed to contain the sum over the rows of the inverse of the sparse matrix. But "out" contains only zeros. Also isstat from MKL_SPARSE_D_TRSM does not indicate any error.

The compiling and linking was:

ifort --version
ifort (IFORT) 19.0.2.187 20190117

ifort -i8 -warn nounused -warn declarations -O3 -static -align array64byte -mkl=parallel -qopenmp -parallel -c -o OMP_MKLPARA_ifort_4.20.12-arch1-1-ARCH/Test.o Test.f90 -I /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/include/

ifort -i8 -warn nounused -warn declarations -O3 -static -align array64byte -mkl=parallel -qopenmp -parallel -o Test_OMP_MKLPARA_4.20.12-arch1-1-ARCH OMP_MKLPARA_ifort_4.20.12-arch1-1-ARCH/Test.o /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/lib/intel64/libmkl_blas95_ilp64.a /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/lib/intel64/libmkl_lapack95_ilp64.a -Wl,--start-group /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/lib/intel64/libmkl_intel_ilp64.a /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/lib/intel64/libmkl_core.a /opt/intel/compilers_and_libraries_2019.2.187/linux/mkl/lib/intel64/libmkl_intel_thread.a -Wl,--end-group -lpthread -lm -ldl

 

Any suggestion?

0 Kudos
7 Replies
may_ka
Beginner
461 Views

I could help myself ......................... trsm is only supported for "SPARSE_MATRIX_TYPE_TRIANGULAR". That means that I have to come up with a popper decomposition of A. Is there any MKL routine which provides that in a format reusable for Inspector-Executor routines??

Thanks

0 Kudos
Gennady_F_Intel
Moderator
461 Views

i am not sure understand your question but may be mkl_sparse_?_trsm may fits with your needs.

0 Kudos
Kirill_V_Intel
Employee
461 Views

Hi,

As the documentation says (and you figured out by yourself), mkl_sparse_?_trsm (and trsv) is designed to work with triangular matrices only. For now we don't have a routine for solving systems with general sparse matrices within Inspector-Executor Sparse BLAS. What you can do is using Intel MKL PARDISO for that and wrap it up for IE Sparse BLAS.

Best,
Kirill

0 Kudos
Kirill_V_Intel
Employee
461 Views

Hi,

A quick follow-up. IE Sparse BLAS is mostly beneficial when used with mkl_sparse_optimize. It might be a good idea to consider adding this optimization step to your code. Of course, unless you have a good reason not to do it.

Best,
Kirill

0 Kudos
may_ka
Beginner
461 Views

Hi, thanks for the elaborate response. I wanted to solve a system y=Ax, where A is sparse and symmetric. A pointed out by Kirill, one can do that with pardiso .......................... and I have done this so far. However, as discussed here there is a problem with pardiso when called repeatedly when the number of right hand sides is small.

I read about trsm when starting to change the code around pardiso to such that I can replace pardiso by a nvidia cuda solver. I must admit that I was not happy to move that way but pardiso is part of my large scale preconditioned gradient solver used to solve up to 500 Mio equations and has become a major bottleneck more than doubling the seconds per iteration. So I would strongly suggest intel to separately supply functions to factorize and subsequently foreward-backward solve matrices (similar to cholmod).

Anyway, with regard to trsm I indeed picked up the problem of feeding a general matrix into trsm after I additionally put in the optimizer routines. These routines reported an error "matrix not supported". however, since trsm can be run without an optimizer call it should also pick up a non-elidgible matrix type and report that.

cheers

0 Kudos
Kirill_V_Intel
Employee
461 Views

Hi,

As Gennady suggested, you can submit a request via https://supporttickets.intel.com/?lang=en-US to speed up performance improvements.

As for the trsm behavior, thanks for the catch!

We'll surely add the format check to this functionality.

Best,
Kirill

 

0 Kudos
Gennady_F_Intel
Moderator
461 Views

note - we added format checking for TRSM. Now, if a matrix with type SPARSE_MATRIX_TYPE_GENERAL is provided as input, the routine will return status SPARSE_STATUS_NOT_SUPPORTED. The feature available since MKL 2019u5 ( which already released the last week).

0 Kudos
Reply