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

knhkse — Thu, 17 Jul 2025 03:42:03 GMT

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

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.

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

noffermans — Fri, 25 Jul 2025 09:38:33 GMT

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

topic Triangular matrix for BSR format - Only supports continuous row pointer array? in Intel® oneAPI Math Kernel Library

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

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