Hi,

Petros · ‎07-16-2015

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!

Ying_H_Intel · ‎07-16-2015

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;
}

Anjan_Mukherjee · ‎02-21-2020

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

mecej4 · ‎02-25-2020

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.

Kirill_V_Intel · ‎02-25-2020

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

Holod__Ihor · ‎03-10-2020

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

yang__dong · ‎04-13-2020

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

Kirill_V_Intel · ‎04-13-2020

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

Holod__Ihor · ‎05-05-2020

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.

Gennady_F_Intel · ‎05-06-2020

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

mkl_?csrcoo duplicate COO values