COO (with duplicates) to CSR

elmeliegy__abdelrahm · ‎05-09-2018

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.

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

mecej4 · ‎05-10-2018

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(N_nz²) operations. When N_nz, 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(N_nz lg N_nz) operations.

What I suggest at this point is that you use MKL_?CSRCOO with job[0] = 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(N_nz) 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.

Alexander_K_Intel2 · ‎05-10-2018

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

mecej4 · ‎05-10-2018

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.