Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- mkl_sparse_sypr value update in array B

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

may_ka

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-05-2020
07:36 AM

381 Views

mkl_sparse_sypr value update in array B

Hi intel mkl team,

from trial and error I inferred that the inspector executor routines "MKL_SPARSE_?_CREATE_?" sets pointers to the supplied arrays instead of taking copies. To my knowledge this is not explicitly mentioned in the manual.

Therefore changing values in the original array directly (NOT via the update function) will affect related operation results (e.g. matrix multiplication). This holds for "mkl_sparse_sp2" and "mkl_sparse_?_mm".

However, this does not seem to hold for "mkl_sparse_?_sypr" where updating values in array B has no effect on the result of "C=A'BA". It appears as if a copy of B is generated when calling "mkl_sparse_?_sypr". A workaround is the destruction of rebuilding of B, with the subsequent requirement to destruct and rebuild C. In some applications, where "C=A'BA" is performed repeatedly and only values in B change but not the non-zero pattern, destroying B and C is an unnecessary overhead.

The code below implements "mkl_sparse_?_sypr" and shows that updating values in B directly has no effect on "C=A'BA". As lined out above this behavior appears inconsistent with other inspector executor multiplication routines and I am wondering whether that is a bug.

Further, it appears as if "mkl_sparse_?_sypr" is not supporting symmetric matrices of storage type "sparse_matrix_type_symmetric". If the code below is run on array "sbfu" the multiplication call throws "SPARSE_STATUS_NOT_SUPPORTED".

thanks for looking into this.

```
#include "mkl_spblas.f90"
module mod_xxx
use mkl_spblas
use, intrinsic :: iso_fortran_env, only: int64, real64,int32
use iso_c_binding, only: c_ptr, c_int, c_long_long, c_double,&
& c_f_pointer
type :: sparse
type(sparse_matrix_t), allocatable :: handle
type(matrix_descr), allocatable :: d
integer(int64), allocatable :: rp(:),cp(:)
real(real64), allocatable :: v(:)
integer(int64), pointer, contiguous :: rpp1(:)=>null(),rpp2(:)=>null(),cpp(:)=>null()
real(real64), pointer, contiguous :: vp(:)=>null()
integer:: nrows, ncols
contains
procedure :: iii => subiii
procedure :: setpointer => subsetpointer
end type sparse
contains
subroutine subiii(this,stat)
implicit none
class(sparse), intent(inout) :: this
integer, intent(inout) :: stat
stat=MKL_SPARSE_D_CREATE_CSR(&
&this%handle,&
&sparse_index_base_one,&
&this%nrows,&
&this%ncols,&
&this%rp(1:(size(this%rp)-1)),&
&this%rp(2:size(this%rp)),&
&this%cp,&
&this%v)
end subroutine subiii
subroutine subsetpointer(this,stat)
implicit none
class(sparse), intent(inout) :: this
integer, intent(inout) :: stat
integer(c_int) :: index
type(c_ptr) :: rows_start,rows_end,col_indx,values
index=sparse_index_base_one
stat=MKL_SPARSE_D_EXPORT_CSR(&
&this%handle,&
&index,&
&this%nrows,&
&this%ncols,&
&rows_start,&
&rows_end,&
&col_indx,&
&values)
if(stat/=0) return
call c_f_pointer(rows_start,this%rpp1,[this%nrows])
call c_f_pointer(rows_end,this%rpp2,[this%nrows])
call c_f_pointer(col_indx,this%cpp,[this%rpp2(this%nrows)-1])
call c_f_pointer(values,this%vp,[this%rpp2(this%nrows)-1])
end subroutine subsetpointer
function err(stat) result(x)
implicit none
integer, intent(in) :: stat
character(:), allocatable :: x
select case(stat)
case(SPARSE_STATUS_NOT_INITIALIZED)
x="SPARSE_STATUS_NOT_INITIALIZED"
case(SPARSE_STATUS_ALLOC_FAILED)
x="SPARSE_STATUS_ALLOC_FAILED"
case(SPARSE_STATUS_INVALID_VALUE)
x="SPARSE_STATUS_INVALID_VALUE"
case(SPARSE_STATUS_EXECUTION_FAILED)
x="SPARSE_STATUS_EXECUTION_FAILED"
case(SPARSE_STATUS_INTERNAL_ERROR)
x="SPARSE_STATUS_INTERNAL_ERROR"
case(SPARSE_STATUS_NOT_SUPPORTED)
x="SPARSE_STATUS_NOT_SUPPORTED"
case default
x="unknown"
end select
end function err
end module mod_xxx
program test
use mkl_spblas
use mod_xxx, only: sparse, err
use, intrinsic :: iso_fortran_env, only: int64, real64,int32
implicit none
type(sparse), allocatable :: sA,sbfu,sC,sbut
real(real64), allocatable :: rA(:,:),rB(:,:),rC(:,:)
integer :: stat=0, i
!!@@@@@@@@@@@@@@@@@
!!the matrix B
!!0.5,0.0,0.1
!!0.0,0.4,0.3
!!0.1,0.3,0.2
!!@@@@@@@@
!!matrix B dense storage
allocate(rA(3,3),rB(3,3),rC(3,3))
rB(1,:)=(/0.5,0.0,0.1/)
rB(2,:)=(/0.0,0.4,0.3/)
rB(3,:)=(/0.1,0.3,0.2/)
!!@@@@@@@@@@
!!the matrix A
!!0.0,1.0,0.0
!!0.0,0.0,1.0
!!0.0,1.0,0.0
!!@@@@@@@@@@@
!!matrix A dense storage
rA(1,:)=(/0.0,1.0,0.0/)
rA(2,:)=(/0.0,0.0,1.0/)
rA(3,:)=(/0.0,1.0,0.0/)
!!multiplication result from dense matrices
rC=matmul(matmul(transpose(rA),rB),rA)
do i=1,size(rC,1)
write(*,"(*(g0:"",""))") rC(i,:)
end do
!!@@@@@@@@@@@@@@@@@@@@@@
allocate(sA,sbfu,sbut,sC);allocate(sA%handle,sbfu%handle,sbut%handle,sC%handle)
allocate(sA%d,sbfu%d,sbut%d,sC%d)
!!@@@@@@@@@@
!!matrix A sparse full storage
sA%d%type=sparse_matrix_type_general
sA%d%mode=SPARSE_FILL_MODE_FULL
sA%d%diag=sparse_diag_non_unit
sA%rp=(/1,2,3,4/)
sA%cp=(/2,3,2/)
sA%v=(/1.0,1.0,1.0/)
sa%nrows=3;sa%ncols=3
call sa%iii(stat)
if(stat/=0) call se("err_1"//err(stat))
!!@@@@@@@@@@@@@@
!!matrix B sparse full storage
sbfu%d%type=sparse_matrix_type_symmetric
sbfu%d%mode=SPARSE_FILL_MODE_FULL
sbfu%d%diag=sparse_diag_non_unit
sbfu%rp=(/1,3,5,8/)
sbfu%cp=(/1,3,2,3,1,2,3/)
sbfu%v=(/0.5,0.1,0.4,0.3,0.1,0.3,0.2/)
sbfu%nrows=3;sbfu%ncols=3
call sbfu%iii(stat)
if(stat/=0) call se("err_2"//err(stat))
!!@@@@@@@@@@@@@@
!!matrix B sparse upper triangular storage
sbut%d%type=sparse_matrix_type_symmetric
sbut%d%mode=SPARSE_FILL_MODE_UPPER
sbut%d%diag=sparse_diag_non_unit
sbut%rp=(/1,3,5,6/)
sbut%cp=(/1,3,2,3,3/)
sbut%v=(/0.5,0.1,0.4,0.3,0.2/)
sbut%nrows=3;sbut%ncols=3
call sbut%iii(stat)
if(stat/=0) call se("err_3"//err(stat))
!!@@@@@@@@@@@@@@@@@@@@@@@@@
!!mkl_sparse_sypr 1. round
stat=mkl_sparse_sypr(&
&operation=sparse_operation_transpose,&
&A=sa%handle,&
&descrB=sbut%d,&
&B=sbut%handle,&
&req=SPARSE_STAGE_FULL_MULT,&
&C=sc%handle)
if(stat/=0) call se("err_4 "//err(stat))
!!set pointers to internal structure of c
call sc%setpointer(stat)
if(stat/=0) call se("err_5 "//err(stat))
write(*,"(*(g0:"",""))") sc%rpp1
write(*,"(*(g0:"",""))") sc%rpp2
write(*,"(*(g0:"",""))") sc%cpp
write(*,"(*(g0:"",""))") sc%vp
!!change values in B but not structure of B
sbut%v(:)=sbut%v(:)*5._real64
!!set pointers to internal structure of b
call sbut%setpointer(stat)
!!values in b have changed
write(*,"(*(g0:"",""))") sbut%vp
!!mkl_sparse_sypr 2. round
stat=mkl_sparse_sypr(&
&operation=sparse_operation_transpose,&
&A=sa%handle,&
&descrB=sbut%d,&
&B=sbut%handle,&
&req=SPARSE_STAGE_FULL_MULT,&
&C=sc%handle)
if(stat/=0) call se("err_6 "//err(stat))
!!results are not different from first round
write(*,"(*(g0:"",""))") sc%vp
!!workaround may require destruction of B
contains
subroutine se(c0e)
character(*), intent(in) :: c0e
write(*,*) c0e
call exit(1)
end subroutine se
end program test
```

Link Copied

5 Replies

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-09-2020
04:09 PM

354 Views

Hi @may_ka,

Let me try to address your questions.

1) "from trial and error I inferred that the inspector executor routines "MKL_SPARSE_?_CREATE_?" sets pointers to the supplied arrays instead of taking copies. To my knowledge this is not explicitly mentioned in the manual."

Right. I agree we should say it explicitly in the documentation.

2) "Therefore changing values in the original array directly (NOT via the update function) will affect related operation results (e.g. matrix multiplication). This holds for "mkl_sparse_sp2m" and "mkl_sparse_?_mm".

However, this does not seem to hold for "mkl_sparse_?_sypr" where updating values in array B has no effect on the result of "C=A'BA". It appears as if a copy of B is generated when calling "mkl_sparse_?_sypr". A workaround is the destruction of rebuilding of B, with the subsequent requirement to destruct and rebuild C. In some applications, where "C=A'BA" is performed repeatedly and only values in B change but not the non-zero pattern, destroying B and C is an unnecessary overhead.

The code below implements "mkl_sparse_?_sypr" and shows that updating values in B directly has no effect on "C=A'BA". As lined out"

Correct and this behavior is expected. Unlike sp2m and mm, the sypr routine creates an internal representation and uses it and hence the changes done for the original input arrays do not affect the result of the second sypr operation. As we recognize a need for this, we added mkl_sparse_?_update_values (https://software.intel.com/content/www/us/en/develop/documentation/mkl-developer-reference-c/top/bla...). I believe currently it might not interact with sypr properly but we can open a FR to support your case.

3) "Further, it appears as if "mkl_sparse_?_sypr" is not supporting symmetric matrices of storage type "sparse_matrix_type_symmetric". If the code below is run on array "sbfu" the multiplication call throws "SPARSE_STATUS_NOT_SUPPORTED"."

Not fully correct, the routine mkl_sparse_?_sypr does not support a combination of sparse_matrix_type_symmetric + SPARSE_FILL_MODE_FULL. It supports symmetric matrices only if they are provided with SPARSE_FILL_MODE_UPPER/LOWER. The documentation says that (in somewhat obscure form but still).

My question here is: do you want the result of the sypr operation also as a full matrix (not just the symmetric part of it)?

Hope this helps.

Best,

Kirill

may_ka

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-11-2020
06:20 PM

340 Views

thanks for the clarification.

With regard to "trial and error" I think till version 17 or 18, inspector-executor routines were taking copies, and with 19 it changed to setting pointers. So an explicit clarification in the manual would be much appreciated as it is pretty important for memory management responsibilities.

I would much appreciate if the "_update" routine would properly interact with with "_sypr". But updating means copying and therefore not taking an internal copy would be even better. However, for performance reason the latter might not be possible. As mentioned C=ABA' is an operation repeated several hundred thousand times, and any unnecessary memory reallocation is an avoidable overhead.

Are you saying that one should feed a full store matrix into the algorithm pretending it is upper/lower stored? If you are NOT saying that allowing only upper/lower storage means building an upper/lower array from a full store array for exactly this operation (assuming that other steps in the process need a full storage). Since some mkl facilities (like pardiso) and other libraries require specific storage formats, it would be good if "_sypr" would support all combinations of upper/lower and full for input and output (e.g. input=lower/upper, output=full etc.). Then the user can optimize the input and output formats with respect to the work flow requirements without unnecessary copying, creation and destruction.

Thank and best regards

Karl

Kirill_V_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-18-2020
07:00 PM

307 Views

Hi @may_ka,

I'm saying that currently mkl_sparse_sypr does not support the case when descrB.type = SPARSE_MATRIX_TYPE_SYMMETRIC and descrB.mode = SPARSE_FILL_MODE_FULL. Only descrB.mode = SPARSE_FILL_MODE_LOWER and SPARSE_FILL_MODE_UPPER are supported (and the documentation says that).

Let us know ( = confirm) if what you need is that mkl_sparse_sypr works with SPARSE_FILL_MODE_FULL.

@Gennady_F_Intel, could you open a FR for extending support of mkl_sparse_update_values to interact nicely with mkl_sparse_sypr?

Best,

Kirill

Steveni

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-11-2020
10:02 PM

335 Views

Hiii,

Thank you for posting in Intel mkl sparse_sypr value update in array B. It is really helpful.

Gennady_F_Intel

Moderator

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

10-19-2020
01:06 AM

299 Views

Ok. We will add this FR to our future plans.

Topic Options

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

For more complete information about compiler optimizations, see our Optimization Notice.