Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
101 Views

mkl_?csrcoo duplicate COO values

Hi,

I am using an old sparse linear solver in my code (of HSL) that receives the matrix in COO form and it accepts to have duplicate values. For example, (1,2,50) and (1,2,30) are two different entries and are summed up internally.

I want to switch to MKL PARDISO to see about the potential speedup I can have from this. From what I get, PARDISO needs the sparse matrix to be in CSR format. Thus, I thought of using mkl_dcsrcoo with job(1)=1 to get my matrix. The problem is, I don't know if it takes care of the duplicates internally or not.

Thanks in advance!

0 Kudos
9 Replies
Highlighted
Employee
101 Views

Hi, 

I guess, the answer is no.   The function can't  handle such sum.  MKL provide example code dconverters.c or f  under MKL example folder. You may try it out. 

and pardiso require the column index are in increasing order, so you may need job(1)=2.

Best Regards,

Ying 

/*

********************************************************************************
*   Copyright(C) 2004-2014 Intel Corporation. All Rights Reserved.
*   
*   The source code, information  and  material ("Material") contained herein is
*   owned  by Intel Corporation or its suppliers or licensors, and title to such
*   Material remains  with Intel Corporation  or its suppliers or licensors. The
*   Material  contains proprietary information  of  Intel or  its  suppliers and
*   licensors. The  Material is protected by worldwide copyright laws and treaty
*   provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
*   modified, published, uploaded, posted, transmitted, distributed or disclosed
*   in any way  without Intel's  prior  express written  permission. No  license
*   under  any patent, copyright  or  other intellectual property rights  in the
*   Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
*   implication, inducement,  estoppel or  otherwise.  Any  license  under  such
*   intellectual  property  rights must  be express  and  approved  by  Intel in
*   writing.
*   
*   *Third Party trademarks are the property of their respective owners.
*   
*   Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
*   this  notice or  any other notice embedded  in Materials by Intel or Intel's
*   suppliers or licensors in any way.

********************************************************************************
*   Content : MKL Sparse format converters C example
*
********************************************************************************
*
* Example program for using MKL Sparse format converters
* The following Sparse  Sparse format converters are used in the example:
*
*          MKL_DDNSCSR
*          MKL_DCSRCOO
*          MKL_DCSRBSR
*          MKL_DCSRDIA
*          MKL_DCSRCSC
*          MKL_DCSRSKY
*
*/
#include <stdio.h>
#include "mkl_spblas.h"
#include "mkl_types.h"


int main()
   {
/********************************************************************************
*     Definition arrays for sparse matrix formats
********************************************************************************/
#define M      4
#define N      4
#define LDA    4
#define NZMAX  8
#define NNZ    8
#define MBLK   2
#define NN     2
#define INFO   0
#define MN     16
#define IBASE1 1
#define IBASE2 1
#define LOCAT  2
#define IDIAG 3
#define NDIAG 4
#define INDIA 12
        MKL_INT        m = M, n=N, lda=LDA, nzmax=NZMAX, nnz = NNZ, mblk=MBLK, nn=NN, info=INFO,mn=MN;
        MKL_INT            ibase1 = IBASE1,ibase2 = IBASE2, locat = LOCAT,idiag = IDIAG,ndiag = NDIAG;
        double        Adns[MN];
        double        Adns_standard[MN];
        double        Absr[NZMAX];
        double        Absr_standard[NZMAX]      =  {5.0, 9.0, 8.0, 2.0, 3.0, 1.0, 6.0, 4.0};
        double        Acsr[NZMAX];
        double        Acsr_standard[NZMAX]      =  {5.0, 8.0, 9.0, 2.0, 3.0, 6.0, 1.0, 4.0};
        double           Acsc[NZMAX];
        double           Acsc_standard[NZMAX]      =  {5.0, 9.0, 8.0, 2.0, 3.0,
                                                  1.0, 6.0, 4.0};
        double           Adia[INDIA];
        double           Adia_standard[INDIA]      =  {0.0, 9.0, 0.0, 1.0,
                                     5.0, 2.0, 3.0, 4.0,
                                    8.0, 0.0, 6.0, 0.0};
        double           Asky[6];
        double           Askyl_standard[6]          =  {5.0, 9.0, 2.0, 3.0, 1.0, 4.0};
        double        Acoo[NZMAX];
        double        Acoo_standard[NZMAX]      =  {5.0, 8.0, 9.0, 2.0, 3.0, 6.0, 1.0, 4.0};
        MKL_INT        AI[M+1];
        MKL_INT        AI1[M+1];
        MKL_INT        AI_standard[M+1]          =  {1, 3, 5, 7, 9};
        MKL_INT        AJ[NZMAX];
        MKL_INT        AJ1[NZMAX];
        MKL_INT        AJ_standard[NZMAX]        =  {1, 2, 1, 2, 3, 4, 3, 4};
        MKL_INT        AJB[NN];
        MKL_INT        AJB_standard[MBLK]        =  {1,2};
        MKL_INT        AIB[NN+1];
        MKL_INT        AIB_standard[MBLK+1]      =  {1, 2, 3};
        MKL_INT        ir[NZMAX];
        MKL_INT        ir_standard[NZMAX]        =  {1, 1, 2, 2, 3, 3, 4, 4};
        MKL_INT        jc[NZMAX];
        MKL_INT        jc_standard[NZMAX]        =  {1, 2, 1, 2, 3, 4, 3, 4};
        MKL_INT        pointers[M+1];
        MKL_INT        pointersl_standard[M+1]    =  {1, 2, 4, 5, 7};
        MKL_INT        distance[IDIAG];
        MKL_INT        distance_standard[IDIAG]  =  {-1, 0, 1};
        MKL_INT        AJL[6];
        MKL_INT        AJL_standard[6]           =  {1, 1, 2, 3, 3, 4};
//*************************************************************************************************
//*    Declaration of local variables :
//*************************************************************************************************
      MKL_INT job0=1;
      MKL_INT ifail=1;
      MKL_INT nr,ldAbsr;
      MKL_INT i,j,ij;
      double  rfail=0.0;
      MKL_INT job[8];
      printf("\n EXAMPLE PROGRAM FOR CONVERTER FROM ONE\n");
      printf("\n    SPARSE FORMAT ROUTINES TO OTHER    \n");
      printf("\n       REAL SINGLE PRECISION           \n");
      locat=2;
      ibase1=1;
      ibase2=1;
      job[1]=ibase1;
      job[2]=ibase2;
      job[3]=locat;
      job[4]=nzmax;
//***************************************************************************************************
//* TASK 1    Obtain compressed sparse row matrix from dense matrix
//**************************************************************************************************
      for ( j=0; j<n; j++)
         for ( i=0; i<m; i++)
               Adns[i + lda*j]=0.0;

      Adns[0]=5.0;
      Adns[1]=9.0;
      Adns[4]=8.0;
      Adns[5]=2.0;
      Adns[10]=3.0;
      Adns[11]=1.0;
      Adns[14]=6.0;
      Adns[15]=4.0;


      for ( j=0; j<n; j++)
         for ( i=0; i<m; i++)
           Adns_standard[i + lda*j]=Adns[i + lda*j];
      job[0]=0;
      job[5]=1;
      mkl_ddnscsr(job,&m,&n,Adns,&lda,Acsr,AJ,AI,&info);
      if (info!=0) goto FAILURE1;
      for ( i=0; i<n+1; i++)
      {
           ifail=AI-AI_standard;
           if (ifail!=0) goto FAILURE1;
      }
      for ( i=0; i<nzmax; i++)
      {
           ifail=AJ-AJ_standard;
           if (ifail!=0) goto FAILURE1;
      }
      for ( i=0; i<nzmax; i++)
      {
           rfail=Acsr-Acsr_standard;
           if (rfail!=0) goto FAILURE1;
      }
//***************************************************************************************************
//* TASK 2    Obtain dense matrix from compressed sparse row matrix
//***************************************************************************************************
      for ( j=0; j<n; j++)
         for ( i=0; i<m; i++)
               Adns[i+j*lda]=0.0;

      job[0]=1;
      mkl_ddnscsr(job,&m,&n,Adns,&lda,Acsr,AJ,AI,&info);
      if (info!=0) goto FAILURE2;
      for ( j=0; j<n; j++){
         for ( i=0; i<m; i++){
               rfail = Adns_standard[i + lda*j]-Adns[i + lda*j];
           if (rfail!=0) goto FAILURE2;
          }
      }
//********************************************************************************************
//* TASK 3    Obtain sparse coordinate matrix from cpompressed sparse row matrix 
//***************************************************************************************************
      job[0]=0;
      job[5]=3;
      mkl_dcsrcoo (job,&n, Acsr, AJ,AI,&nnz,Acoo, ir,jc,&info);
      if (info!=0) goto FAILURE3;
      for ( i=0; i<nzmax; i++)
      {
           ifail=ir-ir_standard;
           if (ifail!=0) goto FAILURE3;
      }
      for ( i=0; i<nzmax; i++)
      {
               ifail=jc-jc_standard;
               if (ifail!=0) goto FAILURE3;
      }

      for ( i=0; i<nzmax; i++)
      {
               rfail=Acoo-Acoo_standard;
           if (rfail!=0) goto FAILURE3;
      }
//***************************************************************************************************
//* TASK 4    Obtain compressed sparse row matrix from sparse coordinate matrix 
//***************************************************************************************************
      job[0]=1;
      job[5]=2;
      mkl_dcsrcoo (job,&n, Acsr, AJ,AI,&nnz,Acoo, ir,jc,&info);
      if (info!=0) goto FAILURE4;
      for ( i=0; i<n+1; i++)
      {
               ifail=AI-AI_standard;
               if (ifail!=0) goto FAILURE4;
      }
      for ( i=0; i<nzmax; i++)
      {
               ifail=AJ-AJ_standard;
               if (ifail!=0) goto FAILURE4;
      }
      for ( i=0; i<nzmax; i++)
      {
               rfail=Acsr-Acsr_standard;
               if (rfail!=0) goto FAILURE4;
      }
//***************************************************************************************************
//* TASK 5    Obtain block sparse row matrix from compressed sparse row matrix 
//***************************************************************************************************
      ldAbsr=mblk*mblk;
      for ( i=0; i < nzmax; i++)
           Absr=0.0;
      job[0] = 0;
      job[2] = 1;
      job[5] = 1;
      mkl_dcsrbsr (job,&m,&mblk,&ldAbsr,Acsr,AJ,AI,Absr,AJB,AIB,&info);

      nr = 1 + (m-1) / mblk;
      if (info!=0) goto FAILURE5;
      for ( i=0; i < nr+1; i++)
      {
               ifail=AIB-AIB_standard;
               if (ifail!=0) goto FAILURE5;
      }
      for ( i=0; i < mblk; i++)
      {
               ifail=AJB-AJB_standard;
               if (ifail!=0) goto FAILURE5;
      }

      for ( i=0; i < nzmax; i++)
      {
               rfail=Absr-Absr_standard;
               if (rfail!=0) goto FAILURE5;
      }
//***************************************************************************************************
//* TASK 6    Obtain compressed sparse row matrix from block sparse row matrix
//***************************************************************************************************
      ldAbsr=mblk*mblk;
      for ( i=0; i<nzmax; i++)
               Acsr=0;
      job[0] = 1;
      job[5] = 3;
      mkl_dcsrbsr (job,&nn,&mblk,&ldAbsr,Acsr,AJ,AI,Absr,AJB,AIB,&info);

      if (info!=0) goto FAILURE6;
      for ( i=0; i< n+1; i++)
      {
               ifail=AI-AI_standard;
               if (ifail!=0) goto FAILURE6;
      }
      for ( i=0; i < nzmax; i++)
      {
              ifail=AJ-AJ_standard;
              if (ifail!=0) goto FAILURE6;
      }
      for ( i=0; i < nzmax; i++)
      {
              rfail=Acsr-Acsr_standard;
              if (rfail!=0) goto FAILURE6;
      }
//***************************************************************************************************
//* TASK 7   Obtain diagonal matrix from compressed sparse row matrix
//***************************************************************************************************
      for ( i=0; i< idiag; i++)
               distance=distance_standard;
      job[0] = 0;
      job[5] = 0;
      mkl_dcsrdia(job,&n,Acsr,AJ,AI,Adia,&ndiag,distance,&idiag,
                  Acsr,AJ,AI,&info);
      if (info!=0) goto FAILURE7;
      for ( i=0; i< idiag; i++)
      {
               ifail=distance-distance_standard;
               if (ifail!=0) goto FAILURE7;
      }
      for ( j=0; j< idiag; j++)
         for ( i=0; i< ndiag; i++)
      {
               ij = i + ndiag*j;
               rfail=Adia-Adia_standard;
               if (rfail!=0) goto FAILURE7;
      }
//***************************************************************************************************
//* TASK 8    Obtain compressed sparse row matrix from  diagonal matrix
//***************************************************************************************************
      job[0] = 1;
      job[5] = 0;

      mkl_dcsrdia(job,&n,Acsr,AJ,AI,Adia,&ndiag,distance,&idiag,Acsr,AJ,AI,&info);

      if (info!=0) goto FAILURE8;
      for ( i=0; i< n+1; i++)
      {
               ifail=AI-AI_standard;
               if (ifail!=0) goto FAILURE8;
      }
      for ( i=0; i < nzmax; i++)
      {
               ifail=AJ-AJ_standard;
               if (ifail!=0) goto FAILURE8;
      }

      for ( i=0; i < nzmax; i++)
      {
               rfail=Acsr-Acsr_standard;
               if (rfail!=0) goto FAILURE6;
      };
//***************************************************************************************************
//* TASK 9    Obtain compressed sparse column matrix from compressed sparse row matrix 
//***************************************************************************************************
      job[0] = 0;
      job[5] = 1;

      mkl_dcsrcsc(job,&n,Acsr,AJ,AI,Acsc,AJ1,AI1,&info);

      if (info!=0) goto FAILURE9;
      for ( i=0; i< n+1; i++)
      {
               ifail=AI1-AI_standard;
               if (ifail!=0) goto FAILURE9;
      }
      for ( i=0; i < nzmax; i++)
      {
               ifail=AJ1-AJ_standard;
               if (ifail!=0) goto FAILURE9;
      }
       for ( i=0; i < nzmax; i++)
      {
               rfail=Acsc-Acsc_standard;
               if (rfail!=0) goto FAILURE9;
      };
//***************************************************************************************************
//* TASK 10    Obtain compressed sparse row matrix from compressed sparse column matrix
//***************************************************************************************************
      job[0] = 1;
      job[5] = 1;
      mkl_dcsrcsc(job,&n,Acsr,AJ,AI,Acsc,AJ1,AI1,&info);

      if (info!=0) goto FAILURE10;
      for ( i=0; i< n+1; i++)
      {
               ifail=AI-AI_standard;
               if (ifail!=0) goto FAILURE10;
      }
      for ( i=0; i < nzmax; i++)
      {
              ifail=AJ-AJ_standard;
              if (ifail!=0) goto FAILURE10;
      }

      for ( i=0; i < nzmax; i++)
      {
               rfail=Acsr-Acsr_standard;
               if (rfail!=0) goto FAILURE10;
      };
//***************************************************************************************************
//* TASK 11    Obtain skyline matrix for lower triangle from compressed sparse row matrix
//***************************************************************************************************
      job[0] = 0;
      job[3] = 0;
      job[4] = 6;
      job[5] = 0;

      mkl_dcsrsky(job,&n,Acsr,AJ,AI,Asky,pointers,&info);

      if (info!=0) goto FAILURE11;
      for ( i=0; i< n+1; i++)
      {
               ifail=pointers-pointersl_standard;
               if (ifail!=0) goto FAILURE11;
      }
      nnz = pointers - pointers[0];
      for ( i=0; i < nnz; i++)
      {
               rfail=Asky-Askyl_standard;
               if (rfail!=0) goto FAILURE11;
      };
//***************************************************************************************************
//* TASK 12   Obtain compressed sparse row matrix for lower triangle from skyline matrix
//***************************************************************************************************
      job[0] = 1;
      job[3] = 0;
      job[5] = 0;

      mkl_dcsrsky(job,&n,Acsr,AJ,AI,Asky,pointers,&info);

      if (info!=0) goto FAILURE12;
      for ( i=0; i< n+1; i++)
      {
               ifail=AI-pointersl_standard;
               if (ifail!=0) goto FAILURE12;
      }
      nnz = pointers - pointers[0];
      for ( i=0; i < nnz; i++)
      {
               ifail=AJ-AJL_standard;
               if (ifail!=0) goto FAILURE12;
      }

      for ( i=0; i < nnz; i++)
      {
              rfail=Acsr-Askyl_standard;
              if (rfail!=0) goto FAILURE12;
      }
//****************************************************************************************************
      printf("\n           All tests passed  \n");
      return 0;
  /* Failure message to print if something went wrong */
FAILURE1: printf("\n Example FAILED to convert from dns to csr...\n");
      return 1;
FAILURE2: printf("\n Example FAILED to convert from csr to dns...\n");
      return 2;
FAILURE3: printf("\n Example FAILED to convert from csr to coo...\n");
      return 3;
FAILURE4: printf("\n Example FAILED to convert from coo to csr...\n");
      return 4;
FAILURE5: printf("\n Example FAILED to convert from csr to bsr...\n");
      return 5;
FAILURE6: printf("\n Example FAILED to convert from bsr to csr...\n");
      return 6;
FAILURE7: printf("\n Example FAILED to convert from csr to dia...\n");
      return 7;
FAILURE8: printf("\n Example FAILED to convert from dia to csr...\n");
      return 8;
FAILURE9: printf("\n Example FAILED to convert from csr to csc...\n");
      return 9;
FAILURE10: printf("\n Example FAILED to convert from csc to csr...\n");
      return 10;
FAILURE11: printf("\n Example FAILED to convert from csr to sky...\n");
      return 11;
FAILURE12: printf("\n Example FAILED to convert from sky to csr...\n");
      return 12;
   }

0 Kudos
Highlighted
101 Views

I have prepared on code which will add up all elements of repeated column hope this will help someone .

 

matrixsize=M value of M X M matrix

subroutine condense(matrixsize,elements,rowindex,columnindex,nnz)
     implicit none
     integer, intent(in) :: matrixsize
     integer, intent(inout) :: nnz,rowindex(2e5),columnindex(2e5)
     double precision, intent(inout) :: elements(2e5)
     integer, allocatable :: ia(:),ja(:)
     double precision, allocatable :: a(:)
     integer :: mm,nn,i,j
     !variable allocation
     allocate(ia(matrixsize+1))
     allocate(ja(nnz))
     allocate(a(nnz))
     !Assembly of stiffness matrix
     mm=1
     ia(1)=1
     nn=1
     do i=1,matrixsize !lopp over row
          do j=1,rowindex(i+1)-rowindex(i) !loop over all elements of row
               if (j .EQ. 1) then !Storing the first element
                    a(mm)=elements(rowindex(i)+j-1) !Storing the nonzero element
                    ja(mm)=columnindex(rowindex(i)+j-1)!storing the column index
                    mm=mm+1
                    nn=nn+1
               else
                    !From second element checking for repeated index
                    if (columnindex(rowindex(i)-2+j) .EQ. columnindex(rowindex(i)+j-1)) then
                         mm=mm-1
                         nn=nn-1
                         a(mm)=a(mm)+elements(rowindex(i)+j-1)
                         ja(mm)=columnindex(rowindex(i)+j-1)
                    else
                         a(mm)=elements(rowindex(i)+j-1)
                         ja(mm)=columnindex(rowindex(i)+j-1)
                    end if
                    mm=mm+1
                    nn=nn+1
               end if
          end do
          ia(i+1)=nn
     end do
     !updating the nnz
     nnz=nn-1
     !Updating the sparse matrix
     elements(:)=0.d0
     elements(1:nnz)=a(:)
     !updating the row index
     rowindex(:)=0
     rowindex(1:matrixsize+1)=ia(1:matrixsize+1)
     !updating the column index
     columnindex(:)=0
     columnindex(1:nnz)=ja(:)
end subroutine

0 Kudos
Highlighted
Black Belt
101 Views

Anjan Mukherjee:

In situations where multiple entries with the same (row,column) occur, such as when assembling an stiffness matrix for a frame structure, such entries often occur out of order. Your code appears to assume that multiple contributions are contiguous, and that the items are already sorted by row and column.

0 Kudos
Highlighted
Employee
101 Views

Hello Anjan,

FYI, in IE SpBLAS we have mkl_sparse_?_create_coo and mkl_sparse_convert_csr. The converter allows duplicates in the coo format. Then you can call mkl_sparse_?_export_csr and get a CSR format back from the matrix handle.

Best,
Kirill

0 Kudos
Highlighted
Beginner
101 Views

Voronin, Kirill (Intel) wrote:

Hello Anjan,

FYI, in IE SpBLAS we have mkl_sparse_?_create_coo and mkl_sparse_convert_csr. The converter allows duplicates in the coo format. Then you can call mkl_sparse_?_export_csr and get a CSR format back from the matrix handle.

Best,
Kirill

When was the consolidation of duplicates introduced? I'm asking since I have different outputs after consecutive calls of the following functions with different versions of MKL (mkl/2019.5 and mkl/2018.1.163)

mkl_sparse_d_create_coo
mkl_sparse_convert_csr
mkl_sparse_order
mkl_sparse_d_export_csr

 

0 Kudos
Highlighted
Beginner
101 Views

Holod, Ihor wrote:

Quote:

Voronin, Kirill (Intel) wrote:

 

Hello Anjan,

FYI, in IE SpBLAS we have mkl_sparse_?_create_coo and mkl_sparse_convert_csr. The converter allows duplicates in the coo format. Then you can call mkl_sparse_?_export_csr and get a CSR format back from the matrix handle.

Best,
Kirill

 

 

When was the consolidation of duplicates introduced? I'm asking since I have different outputs after consecutive calls of the following functions with different versions of MKL (mkl/2019.5 and mkl/2018.1.163)

mkl_sparse_d_create_coo
mkl_sparse_convert_csr
mkl_sparse_order
mkl_sparse_d_export_csr

 

Hello,

Now I also encounter the same problem. I am sure that the old version function 'mkl_?csrcoo' doesn't add the duplicate values automatically. But when I try to use the new version function 'mkl_sparse_d_create_coo, mkl_sparse_d_convert_csr, mkl_sparse_d_export_csr', the mkl_sparse_d_export_csr always produce wrong results. 

Could you tell me if you have resolved this problem? which version of MKL is ok now?

 

Thanks in advance.

 

Eric

0 Kudos
Highlighted
Employee
101 Views

Hello Eric,

I don't know the answer the question about when.
 

Could you be more specific about when exactly you see wrong results? I'd like to check what is going on there and if we have a bug in MKL, we definitely want to fix it. It'd be great if you can provide us an example with data.

Best,
Kirill

0 Kudos
Highlighted
Beginner
101 Views

I've also noticed that mkl_sparse_d_create_coo doesn't sum up duplicates in version 2018.2.199. Looks like in 2018.4 the summation is in place.

0 Kudos
Highlighted
Moderator
101 Views

in an ideal case, it would be better if you will give us the reproducer and all needed specific details which we use to investigate the issue behavior on our side with any versions of mkl ( including the latest).

0 Kudos