0

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
Alex Reinking
  • 16,724
  • 5
  • 52
  • 86
  • I would check your compiler's documentation on whether C++ 11's range-based `for` loops work reliably on non-standard arrays. Your `data`, `cols` and `rows` are non-standard, as you cannot declare array sizes using a runtime value in standard C++. – PaulMcKenzie Aug 13 '14 at 01:38
  • I'm using icc 14.0.3. If you count the entries printed by the range-based `for` loops, you'll see that it's printing the correct number of elements in all three cases. Just to be sure, though, I rewrote them as normal `for` loops, and nothing (except the syntax) changed. – Alex Reinking Aug 13 '14 at 01:43

0 Answers0