[C Interface] Updating the values of a CSC matrix

[EnricoU]

 Hi!

for my master's thesis I am implementing a simulation routine, which involves for each iteration:

a) matrix-matrix multiplication of a diagonal matrix with a sparse matrix (mkl_sparse_sp2m)

b) Some custom accumulation results from a) yielding a dense vector

c) sparse matrix-vector multiplication with the vector from b). (mkl_sparse_d_mv)

I would like to minimize overhead in this critical portion of the code. See below the message for the code snippets in question.

After some trial and error I have two questions, which the docs could not resolve for me:

1.

The Diagonal matrix xuDiag from a) stays structurally constant. For every iteration, I just want to update  the main diagonal to the result of c). I was hoping to just overwrite the pointer to the values (given at initial creation to mkl_sparse_d_create_csc). This does not seem to have an effect. Since mkl_sparse_?_update_values only supports bsr, is my only option to recreate the whole matrix for every iteration?

I would very much like to avoid this kind of overhead and just write to the internal value array directly or at least have a way of updating only the values without recreating the whole structure for every iteration.

Edit:

After browsing through more related forum Posts I read Spencer_P_Intel's Post, which touches on this topic I believe. This strengthens my first suspicion, that my naive approach of changing the pointer "below the libraries feet" is ill-concieved.

This does leave me with the question, what possibilties do I have for editing values of matrices?

I rely on CSC (CSR in a pinch) for efficient intermediary calculations, so i can only 

• recreate the diagonal matrix for every iteration
• loop through every nnz one by one via mkl_sparse_d_set_value

?

1.5

After reading through Spencer's above mentioned post, there is another question that is on my mind the whole way through learning my way around this library:

Where would I go and learn about these implicit (or explicit) agreements that are mentioned here? I see that in the C++ docs, there seem to be more hints about this I don't expect this information to translate to C, at least not fully?

The same goes for informations like this. It's very useful to know that CSR is the preferred format to use, but how could I have known this prior to seeing the limitations on the seperate function or this forum?

I want to stress:This is not meant as criticism in any way - on the contrary I a very grateful for the great software and maintenance provided by this team. All the more motivation for me ( and maybe future readers) to find that info ourselves and save time for everyone.

2. 

Reading through the documentation, I noticed that the interfaces for C and C++ show noticable differences. For C the Inspector-Executor matrix-matrix multiplication for example only supports the SPARSE_MATRIX_TYPE_GENERAL, while the C++ documentation indicates no such restriction.

Is there a functionally a difference between the C and C++ interfaces, and would one or the other be more advisable in terms of performance regarding sparse inspector-executor operations?

I am very grateful for any input and hope you have a great day.

EnricoU

If it would be of any help, I would gladly attach the sourcecode.

sparse_status_t xuDiag_status = mkl_sparse_d_create_csc(
        &xuDiag,
        SPARSE_INDEX_BASE_ONE,
        U_m, U_m,       // !! U_m x U_m
        nnz_rows,
        (nnz_rows + 1), 
        cols,       
        xu
    );

sparse_status_t Uxu_status = mkl_sparse_sp2m(
        SPARSE_OPERATION_NON_TRANSPOSE,
        descA,
        xuDiag,
        SPARSE_OPERATION_NON_TRANSPOSE,
        descB,
        U,
        SPARSE_STAGE_FULL_MULT,
        &Uxu
);

 sparse_status_t expStatus = mkl_sparse_d_export_csc(
        Uxu,
        &index_base,
        &Uxu_rows,
        &Uxu_cols,
        &Uxu_rows_start,
        &Uxu_rows_end,
        &Uxu_col_inds,
        &Uxu_vals
);

/* Some calculations on Uxu_vals */

const int t_num = 10; 
double **x_result = (double**) mkl_malloc(t_num * sizeof(double*),64);
for (size_t i = 0; i < 10; i++)
{
    x_result[i] = (double*) mkl_malloc(Phi_m * sizeof(double),64);
}

sparse_status_t finalStatus = mkl_sparse_d_mv(
    SPARSE_OPERATION_NON_TRANSPOSE,
    1.0,
    Phi,
    descA,
    Uxucum,
    0.0,        // 1.0 would include the values already being there
    x_result[0]
);

xu = x_result[0];   // raw pointer tomfoolery
vdPackI(
    URef_n,         // number of elements
    URef_vals + 1,  // source
    URef_m,         // step size source. CSC 
    (xu + XRef_n)   // dest  with offset, since u has to come "after" x
);
/* This operation I would like to avoid, since this portion of the code will the most often */

xuDiag_status = mkl_sparse_d_create_csc(
        &xuDiag,
        SPARSE_INDEX_BASE_ONE,
        U_m, U_m,       // !! U_m x U_m
        nnz_rows,
        (nnz_rows + 1), 
        cols,       // For diagonal matrices, rows_start == rows_end  == cols
        xu
);
/* This operation does not work, without re-creating xuDiag as above
*  Internal representation of xuDiag seems to be unimpressed by the reassignment of xu   
*/
Uxu_status = mkl_sparse_sp2m(
    SPARSE_OPERATION_NON_TRANSPOSE,
    descA,
    xuDiag,
    SPARSE_OPERATION_NON_TRANSPOSE,
    descB,
    U,
    SPARSE_STAGE_FINALIZE_MULT, // Only recompute the values, reuse everything else. 
    &Uxu
);

[Aleksandra_K]

Hi Enrico,

Thank you for posting the issue. Please share the code, and I'll look into it more deeply

Regards,

Alex

[Aleksandra_K]

Hi, do you have any updates regarding the code?

[Aleksandra_K]

With no response from you, this issue will no longer be tracked by Intel.

If you have further queries, please feel free to post them in the forum.