Community
cancel
Showing results for 
Search instead for 
Did you mean: 
may_ka
Beginner
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

 

 

0 Kudos
5 Replies
Kirill_V_Intel
Employee
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
340 Views

Hi @Kirill_V_Intel 

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
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
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
299 Views

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

 

Reply