- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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}^{2}) 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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

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 5 4 5] (the fifth value should be **4****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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page