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) 22.214.171.124 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
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??
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.
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.
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.
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.
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).