Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## COO (with duplicates) to CSR Beginner
518 Views
1. HI, I wrote the following routine to convert the e-vectors (in coo with multiple duplicates) to the csr format. The steps I made are (1-sorting ia 2- sorting ja 3- condensing or summing up the duplicates(that val sharing the same(ia&ja) indices 4- Convert using the mkl_csrcoo function. This works fine for a small test example. however it takes too much time for big example.
2. Can anyone help?
```***HERE IS THE CODE***
void Sort_COO(MKL_INT nnz, std::vector<MKL_INT> &ia, std::vector<MKL_INT> &ja, std::vector<double> &val) {
//this part to sort the arrays according to (row and column) array
//first , sort the row array in ascending order and sort the other two arrays corespondigly
//second, sort the column array such that it is ascending within constant row number(i.e within specific row) and sort the other arrays correspondigly.
MKL_INT  holdia, holdja;
double holdVal;
for (int i = 0; i < nnz; i++) {
for (int j = 0; j < nnz - 1; j++) {
if (ia > ia[j + 1]) {
//swapping elements w.r.t row rearrangement
holdia = (ia);
holdja = (ja);
holdVal = (val)
ia = ia[j + 1];  //check this if it works properly ( those are pointers not values)
ja = ja[j + 1];
val = val[j + 1];
(ia)[j + 1] = holdia;
(ja)[j + 1] = holdja;
(val)[j + 1] = holdVal;
}
}
}
for (int i = 0; i < nnz; i++) {
for (int j = 0; j < nnz - 1; j++) {
if ((ja > ja[j + 1]) && (ia == ia[j + 1])) {
//swapping elemnents w.r.t column rearrangement (in assending oreder)
holdia = (ia);
holdja = (ja);
holdVal = (val);
ia = ia[j + 1];
ja = ja[j + 1];
val = val[j + 1];
(ia)[j + 1] = holdia;
(ja)[j + 1] = holdja;
(val)[j + 1] = holdVal;
}
}
}
}
void Condense_COO(MKL_INT& nnz, std::vector<MKL_INT> &ia, std::vector<MKL_INT> &ja, std::vector<double> &val) {
//now, val,ia and ja are constructed in coo with duplicate entriees
//we need to condense this to coo format (without duplications)
//we will sum those values in val which have the same (ia and ja) indecies
MKL_INT nduplic = 0;
for (int i = 0; i < nnz; i++) {
if (((ia) == (ia)[i - 1]) && ((ja) == (ja)[i - 1]))
{
nduplic++;
(val)[i - 1] = (val) + (val)[i - 1];
for (int k = i; k < nnz; k++) {
if (((ia)[k + 1] == (ia)[i - 1]) && ((ja)[k + 1] == (ja)[i - 1]))
{
/*nduplic++;*/
(val)[i - 1] = (val)[k + 1] + (val)[i - 1];
std::cout << "nduplic=" << nduplic << endl;
for (int j = k + 1; j < nnz; ++j)
{
val = val[j + 1];
ia = ia[j + 1];
ja = ja[j + 1];
}
}
}
for (int l = i; l < nnz; ++l)
{
val = val[l + 1];
ia = ia[l + 1];
ja = ja[l + 1];
}
}
}
nnz = nnz - nduplic;
std::cout << "nduplic=" << nduplic << endl;
//erase the surplus elements at the end of each array
for (int i = 0; i < nduplic; ++i)
{
val.pop_back();
ia.pop_back();
ja.pop_back();
}
}
void Convert_CSR(MKL_INT &job, int& N, double &ASCR, MKL_INT& AII, MKL_INT& AJJ, MKL_INT& NNZ, double& VAL, MKL_INT& IA, MKL_INT& JA, MKL_INT& info) {
mkl_dcsrcoo(&job, &N, &ASCR, &AJJ, &AII, &NNZ, &VAL, &IA, &JA, &info);
}```

3 Replies Black Belt
518 Views

When you post code, please use the {...} button on the toolbar. That will make it easier for others to read your code with the indenting preserved and with some syntax highlighting.

You chose bubble sort as the algorithm to sort the matrix entries, and you sorted twice, once on the row and then on the column indices. In addition, you used a similar algorithm to consolidate the multiple entries after sorting. Thus, you have three double loops that take O(Nnz2) operations. When Nnz, the number of non-zero entries, is small, using an inefficient, but simple-to-implement, algorithm is OK. The trouble with this is that when you scale up your problem the inefficiency bites you. You can use a better sorting algorithm such as quicksort or mergesort to do the work in O(Nnz lg Nnz) operations.

What I suggest at this point is that you use MKL_?CSRCOO with job = 2 to do the conversion from COO to CSR. Those MKL routines will do the sorting for you.

However, the output CSR matrix will contain repeated entries (same row and column, with more than one contribution to the value of the matrix entry). You can now consolidate these entries with an O(Nnz) operation as follows. Keep two pointers P1 and P2 to the old (unconsolidated) and new (consolidated) CSR entries. Run through the list and consolidate. When an entry has the same row and column as the previous entry, increment P1, consolidate into the entry pointed by P2; do not increment P2. When the row/column change, increment P1,  store the entry. and update P2.

The current MKL release labels MKL_?CSRCOO as "deprecated", and provides a set of "inspector-executor" routines. I have not been able to get these to work (in Fortran), and I found their documentation to be not yet adequate. There are some code examples in C, but none in Fortran for converting from COO to CSR and then exporting the CSR matrix. If there is a way of consolidating multiple contributions, I have not found it. Employee
518 Views

Hi,

That's correct, MKL inspector-executor functionality handle correctly elements on same position. I add fortran based example that converts rectangular matrix from COO to CSR format

```        PROGRAM test
USE ISO_C_BINDING
USE MKL_SPBLAS

!Declare new interface with aarrays rows, cols, values of type C pointer
!Change name of the routine to EXPORT_CSR_FIXED so it doesn’t conflict with existing MKL interface
INTERFACE
FUNCTION MKL_SPARSE_S_EXPORT_CSR_F(source,indexing,rows,cols,rows_start,rows_end,col_indx,values) &
BIND(C, name='MKL_SPARSE_S_EXPORT_CSR')
USE, INTRINSIC :: ISO_C_BINDING , ONLY : C_INT, C_FLOAT, C_PTR
IMPORT SPARSE_MATRIX_T
TYPE(SPARSE_MATRIX_T), INTENT(IN) :: source
INTEGER(C_INT)      , INTENT(INOUT) :: indexing
INTEGER(C_INT)      , INTENT(INOUT) :: rows
INTEGER(C_INT)      , INTENT(INOUT) :: cols
TYPE(C_PTR)         , INTENT(INOUT) :: rows_start
TYPE(C_PTR)         , INTENT(INOUT) :: rows_end
TYPE(C_PTR)         , INTENT(INOUT) :: col_indx
TYPE(C_PTR)         , INTENT(INOUT) :: values
INTEGER(C_INT) MKL_SPARSE_S_EXPORT_CSR_FIXED
END FUNCTION
END INTERFACE

type (SPARSE_MATRIX_T) :: A_coo,A_csr
integer :: stat
real :: values_a(10)
integer :: col_a(10),row_a(10)
integer :: rows, cols,indexing

!declare C and fortran pointers. C pointers will be passed to the body of export_csr routine.
integer nnz, nrows
type(c_ptr) :: row_start_a_csr, row_end_a_csr, col_a_csr, values_a_csr

integer, pointer, dimension(:) :: row_start_f, row_end_f, col_f
real, pointer, dimension(:) :: values_f

values_a(1:10) = real((/1,2,3,4,5,6,7,8,9,10/))
col_a(1:10) = (/1,2,2,3,4,3,4,5,4,5/)
row_a(1:10) = (/1,2,2,2,2,3,3,3,4,4/)

stat = mkl_sparse_s_create_coo(A_coo,SPARSE_INDEX_BASE_ONE,4,5,10,&
row_a,col_a,values_a)
print *, stat

stat = mkl_sparse_convert_csr(A_coo, SPARSE_OPERATION_NON_TRANSPOSE, A_csr)

!call export_csr with C pointers
stat = mkl_sparse_s_export_csr(A_csr,indexing,rows,cols,&
row_start_a_csr,row_end_a_csr,col_a_csr,values_a_csr)

print *, stat
print *, 'rows ', rows
print *, 'indexing', indexing

nrows = rows

!Cast C pointers that came from export routine to Fortran pointers
call c_f_pointer(row_start_a_csr, row_start_f, [nrows+1])
call c_f_pointer(row_end_a_csr, row_end_f, [nrows])

nnz = row_start_f(nrows+1)-indexing
print *, nnz
call c_f_pointer(col_a_csr,     col_f,     [nnz])
call c_f_pointer(values_a_csr,   values_f,   [nnz])

print *,'row_start_f', row_start_f
print *,'row_end_f', row_end_f
print *,'col_f', col_f
print *,'values_f', values_f

END PROGRAM test

``` Black Belt
518 Views

The code of #3 should show `stat` `= mkl_sparse_s_export_csr_f `on line 49.

The code can be compiled and run with this change, but the results are wrong. With Parallel Studio 18.0.2 on Windows I get

col_f = [1 2 3 4 4 4 5 4 5] (the fifth value should be 3)

values_f = [1 5 4 0 5 7 8 9 10] should be [1 5 4 5 6 7 8 9 10]

Apart from these errors, which you will probably fix (or show me that I made a mistake), I think that having to declare interfaces to C routines, using pointers, ISO_C_BINDING and creating false names for functions will all be major obstacles to invoking these useful matrix manipulation routines from Fortran. After Fortran 95 arrived, most Fortran users use allocatable arrays instead of using the troublesome pointer arrays. Until these concerns can be addressed satisfactorily, please do not remove the deprecated MKL_?CSRCOO routines.

The documentation page https://software.intel.com/en-us/mkl-developer-reference-fortran-mkl-sparse-convert-csr lists "sparse_status_t" as an "output parameter". The correct name, which can be seen at the top of the page, is "stat", and that is the function return value. Furthermore, the word "parameter" is confusing to Fortran users, since in that language a parameter is what is called a "named constant" in other contexts; "output argument" would be better. 