Triangular matrix for BSR format - Only supports continuous row pointer array?

[knhkse]

Hi!

I have a question about `mkl_sparse_d_trsv` routine.

I usually use BSR format, so I want the triangular solving process using each sub-blocks as a unit.

However, I found that only sparse matrix handle with continuous `pointerB` and `pointerE` arrays gives a right solution.

When I used separated arrays for `pointerB` and `pointerE`, the result gave me `SPARSE_STATUS_NOT_SUPPORTED`.

Sorry for my poor English skill, so here is my simple C code for details.

I compiled with

`icx -L${MKLROOT}/lib -lmkl_rt -lm -ldl sample.c -o sample.x` command.

<sample.c>

#include <stdio.h>
#include "mkl.h"

int main(int argc, char *argv[])
{
    sparse_matrix_t op1, op2;
    sparse_status_t stat;
    struct matrix_descr ldescr;
    int bn, blk;

    bn = 3;     // Block row length
    blk = 2;    // Block size

    /**
     * Simple BSR sparse matrix arrays
     * |1, 0, 6, 7, 0, 0|
     * |2, 1, 8, 2, 0, 0|
     * |0, 0, 1, 4, 0, 0|
     * |0, 0, 5, 1, 0, 0|
     * |0, 0, 4, 3, 7, 2|
     * |0, 0, 0, 8, 0, 5|
     */
    int indptr[4] = {0, 2, 3, 5};       // Continuous array
    int pointerB[3] = {0, 2, 3};        // Separated array
    int pointerE[3] = {2, 3, 5};
    int indices[5] = {0, 1, 1, 1, 2};   // Columns
    double data[20] = {                 // Non-zero values
       1., 0., 2., 1., 6., 7., 8., 2., 1., 4., \
       5., 1., 4., 3., 0., 8., 7., 2., 0., 5.
    };

    // Example arrays
    double x[6] = {1., 2., 3., 4., 5., 6.};
    double y[6] = {0., 0., 0., 0., 0., 0.};

    /**
     * 1. Creates BSR handle
     * `op1` : Use continous array (indptr)
     * `op2` : Use separated arrays (pointerB, pointerE)
     */
    stat = mkl_sparse_d_create_bsr(
        &op1, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, \
        bn, bn, blk, indptr, &indptr[1], indices, data
    );
    if (stat != SPARSE_STATUS_SUCCESS) return -1;

    stat = mkl_sparse_d_create_bsr(
        &op2, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, \
        bn, bn, blk, pointerB, pointerE, indices, data
    );
    if (stat != SPARSE_STATUS_SUCCESS) return -1;

    /**
     * 2. Creates matrix description
     * `ldescr` : Block lower triangular matrix
     */
    ldescr.type = SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR;
    ldescr.mode = SPARSE_FILL_MODE_LOWER;
    ldescr.diag = SPARSE_DIAG_NON_UNIT;

    /**
     * 3. Solves linear system with BSR handles
     * `op1` : stat = 0, expected result
     * `op2` : stat = 6, Error occured
     */
    stat = mkl_sparse_d_trsv(
        SPARSE_OPERATION_NON_TRANSPOSE, 1.0, \
        op1, ldescr, x, y
    );
    printf("Status for `op1` : %d\n", stat);
    for (int i=0; i<6; i++) printf("%.4f ", y[i]);
    printf("\n");

    for (int i=0; i<6; i++) y[i] = 0.0;
    stat = mkl_sparse_d_trsv(
        SPARSE_OPERATION_NON_TRANSPOSE, 1.0, \
        op2, ldescr, x, y
    );
    printf("Status for `op2` : %d\n", stat);
    for (int i=0; i<6; i++) printf("%.4f ", y[i]);
    printf("\n");

    return 0;
}

<Output>

Status for `op1` : 0
1.0000 0.0000 0.6842 0.5789 -0.0030 0.2737
Status for `op2` : 6
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

I didn't test every sparse BLAS routines using this circumstance, but the trsv routine is the obvious problem I encountered.

Is there any reason for this? Or, am I missing something?

Thank you.

[noffermans · Accepted solution]

Hi knhkse,

Thank you for reaching out. You are not missing anything and the combination of a BSR matrix in 4-array representation and of a SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR descriptor type is currently not supported.

Given an input matrix in BSR4 format, one possible workaround is to call the   mkl_sparse_convert_bsr() routine, which will internally convert the matrix to BSR3 format. The modified code below (where op3 is op2 converted to BSR3):

#include <stdio.h>
#include "mkl.h"

int main(int argc, char *argv[])
{
    sparse_matrix_t op1, op2, op3;
    sparse_status_t stat;
    struct matrix_descr ldescr;
    int bn, blk;

    bn = 3;     // Block row length
    blk = 2;    // Block size

    /**
     * Simple BSR sparse matrix arrays
     * |1, 0, 6, 7, 0, 0|
     * |2, 1, 8, 2, 0, 0|
     * |0, 0, 1, 4, 0, 0|
     * |0, 0, 5, 1, 0, 0|
     * |0, 0, 4, 3, 7, 2|
     * |0, 0, 0, 8, 0, 5|
     */
    int indptr[4] = {0, 2, 3, 5};       // Continuous array
    int pointerB[3] = {0, 2, 3};        // Separated array
    int pointerE[3] = {2, 3, 5};
    int indices[5] = {0, 1, 1, 1, 2};   // Columns
    double data[20] = {                 // Non-zero values
       1., 0., 2., 1., 6., 7., 8., 2., 1., 4., \
       5., 1., 4., 3., 0., 8., 7., 2., 0., 5.
    };

    // Example arrays
    double x[6] = {1., 2., 3., 4., 5., 6.};
    double y[6] = {0., 0., 0., 0., 0., 0.};

    /**
     * 1. Creates BSR handle
     * `op1` : Use continous array (indptr)
     * `op2` : Use separated arrays (pointerB, pointerE)
     */
    stat = mkl_sparse_d_create_bsr(
        &op1, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, \
        bn, bn, blk, indptr, &indptr[1], indices, data
    );
    if (stat != SPARSE_STATUS_SUCCESS) return -1;

    stat = mkl_sparse_d_create_bsr(
        &op2, SPARSE_INDEX_BASE_ZERO, SPARSE_LAYOUT_ROW_MAJOR, \
        bn, bn, blk, pointerB, pointerE, indices, data
    );
    if (stat != SPARSE_STATUS_SUCCESS) return -1;

    /**
     * Convert op2 in BSR4 format to op3 in BSR3 format
     */
    stat = mkl_sparse_convert_bsr(op2, blk, SPARSE_LAYOUT_ROW_MAJOR, \
                                  SPARSE_OPERATION_NON_TRANSPOSE, &op3);
    if (stat != SPARSE_STATUS_SUCCESS) return -1;
 /**
     * 2. Creates matrix description
     * `ldescr` : Block lower triangular matrix
     */
    ldescr.type = SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR;
    ldescr.mode = SPARSE_FILL_MODE_LOWER;
    ldescr.diag = SPARSE_DIAG_NON_UNIT;

    /**
     * 3. Solves linear system with BSR handles
     * `op1` : stat = 0, expected result
     * `op2` : stat = 6, Error occured
     * `op3` : stat = 0, expected result
     */
    stat = mkl_sparse_d_trsv(
        SPARSE_OPERATION_NON_TRANSPOSE, 1.0, \
        op1, ldescr, x, y
    );
    
    printf("Status for `op1` : %d\n", stat);
    for (int i=0; i<6; i++) printf("%.4f ", y[i]);
    printf("\n");
    
    for (int i=0; i<6; i++) y[i] = 0.0;
    stat = mkl_sparse_d_trsv(
        SPARSE_OPERATION_NON_TRANSPOSE, 1.0, \
        op2, ldescr, x, y
    );
    printf("Status for `op2` : %d\n", stat);
    for (int i=0; i<6; i++) printf("%.4f ", y[i]);
    printf("\n");

    for (int i=0; i<6; i++) y[i] = 0.0;
    stat = mkl_sparse_d_trsv(
        SPARSE_OPERATION_NON_TRANSPOSE, 1.0, \
        op3, ldescr, x, y
    );
    printf("Status for `op3` : %d\n", stat);
    for (int i=0; i<6; i++) printf("%.4f ", y[i]);
    printf("\n");    

    return 0;
}

will produce the following output:

Status for `op1` : 0
1.0000 0.0000 0.6842 0.5789 -0.0030 0.2737
Status for `op2` : 6
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
Status for `op3` : 0
1.0000 0.0000 0.6842 0.5789 -0.0030 0.2737

Hope this helps.

Best,
Nicolas