I'm trying to use the MKL routine mkl_dcsradd
to add an upper-triangular matrix to its transpose. In this case, the upper triangular matrix stores part of the adjacency matrix of a graph, and I need the full version for implementing another algorithm.
In this simplified example, I start with a list of (11) edges, and build an upper-triangular CSR matrix from it. I have checked that this much works. However, when I try to add it to its transpose, dcsradd
stops on the final row, saying it's run out of space. However, this shouldn't be the case. An upper triangular matrix (no zeros along the diagonal) with n
non-zero entries, when added to its transpose, should result in a matrix with 2n
(22) non-zeros.
When I supply dcsradd
with a maximum non-zeros of 22
, it fails, but when I supply it with 23
(an excessive value), it works correctly. Why is this?
I've simplified my code down to a minimal example demonstrating the error:
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <mkl.h>
int main()
{
int nnz = 11;
int numVertices = 10;
int32_t u[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1 };
int32_t v[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 5, 8 };
double w[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
int fullNnz = nnz * 2;
int dim = numVertices;
double triData[nnz];
int triCols[nnz];
int triRows[dim];
// COO to upper-triangular CSR
int info = -1;
int job [] = { 2, 1, 0, 0, nnz, 0 };
mkl_dcsrcoo(job, &dim,
triData, triCols, triRows,
&nnz, w, u, v,
&info);
printf("info = %d\n", info);
// Allocate final memory
double data[fullNnz];
int cols[fullNnz];
int rows[dim];
// B = A + A^T (to make a full adjacency matrix)
int request = 0, sort = 0;
double beta = 1.0;
int WRONG_NNZ = fullNnz + 1; // What is happening here?
mkl_dcsradd("t", &request, &sort, &dim, &dim,
triData, triCols, triRows,
&beta, triData, triCols, triRows,
data, cols, rows,
&WRONG_NNZ, &info);
printf("info = %d\n", info);
// Convert back to 0-based indexing (via Cilk)
cols[:]--;
rows[:]--;
printf("data:");
for (double d : data) printf("%.0f ", d);
printf("\ncols:");
for (int c : cols) printf("%d ", c);
printf("\nrows:");
for (int r : rows) printf("%d ", r);
printf("\n");
return 0;
}
I compile with:
icc -O3 -std=c++11 -xHost main.cpp -o main -openmp -L/opt/intel/composerxe/mkl/lib -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -lpthread -lm
When I give 22, the output is:
info = 0
info = 10
data:1 10 1 2 11 2 3 3 4 4 5 10 5 6 6 7 7 8 11 8 9 0
cols:1 5 0 2 8 1 3 2 4 3 5 0 4 6 5 7 6 8 1 7 9 -1
rows:0 2 5 7 9 11 14 16 18 21
But, when I give 23, the output is:
info = 0
info = 0
data:1 10 1 2 11 2 3 3 4 4 5 10 5 6 6 7 7 8 11 8 9 9
cols:1 5 0 2 8 1 3 2 4 3 5 0 4 6 5 7 6 8 1 7 9 8
rows:0 2 5 7 9 11 14 16 18 21