topic Hi Tianxiong, Mecej4, in Intel® oneAPI Math Kernel Library

What does the means of struct matrix_descr in mol_sparse_XXX routines?

Tianxiong_Lu — Thu, 17 Sep 2015 07:32:52 GMT

Hi all,

In many mkl_sparse_XXX routines (version 11.3), there is a a parameter of type matrix_descr, just like :

sparse_status_t mkl_sparse_s_trsv (sparse_operation_t operation, float alpha, const sparse_matrix_t A, struct matrix_descr descr, const float *x, float *y);

The explanation of this parameter in manual is:

descr : Structure specifying sparse matrix properties. ......But it is confusing. Following is the code snippet in MKL example sparse_trsv.c,

   //*******************************************************************************
    //     Declaration and initialization of parameters for sparse representation of
    //     the matrix A in the compressed sparse row format:
    //*******************************************************************************
#define M 5
#define N 5
#define NNZ 13
    //*******************************************************************************
    //    Sparse representation of the matrix A
    //*******************************************************************************
    double csrVal[NNZ]    = { 1.0, -1.0,     -3.0,
                             -2.0,  5.0,
                                         4.0, 6.0, 4.0,
                             -4.0,       2.0, 7.0,
                                    8.0,          -5.0 };
    MKL_INT    csrColInd[NNZ] = { 0,      1,        3,
                              0,      1,
                                           2,   3,   4,
                              0,           2,   3,
                                      1,             4 };
    MKL_INT    csrRowPtr[M+1] = { 0, 3, 5, 8, 11, 13 };
    // Descriptor of main sparse matrix properties
    struct matrix_descr descrA;
    // Structure with sparse matrix stored in CSR format
    sparse_matrix_t       csrA;

......

    // Create matrix descriptor
    descrA.type = SPARSE_MATRIX_TYPE_TRIANGULAR;
    descrA.mode = SPARSE_FILL_MODE_LOWER;
    descrA.diag = SPARSE_DIAG_UNIT;

    // Compute y = alpha * A^{-1} * x
    mkl_sparse_d_trsv ( SPARSE_OPERATION_NON_TRANSPOSE,
                        alpha,
                        csrA,
                        descrA,
                        x,
                        y );

Obviously, the matrix csrA is not a triangular matrix, and the diagnoal is not unit. Why set the properties of descrA to TRIANGULAR and DIAG_UNIT?
If I change the properties to other option, for example:

    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;
    descrA.diag = SPARSE_DIAG_NON_UNIT;

Now it can not get correct result.

How to setting the properties of matrix_descr correctly?

Best regards,

Tianxiong Lu

Please read the comments in

mecej4 — Thu, 17 Sep 2015 16:34:37 GMT

Please read the comments in the example file carefully. The examples shows how to extract the lower triangular part L of a square matrix A and do things with L. Note that L is stated to have unit diagonal elements, so only the strictly-lower triangular part of A is used. Here are the relevant comment lines:

! The test performs the operation
!
! L^{-1}*x = y, using mkl_sparse_d_trsv
!
! then checks the result:
! 1. compute:
! L*y = y1, using mkl_sparse_d_mv
!
! 2. compare x with y1 - x should be equal y1 (x = L*L^{-1}*x)
!
! Here A is a general sparse matrix
! L is the lower triangular part of matrix A with unit diagonal
! x, y and y1 are vectors
!

Note that every matrix operation is done with L, rather than with A.

If, in your application, you have a different set of needs, and you do not need to extract the lower triangular part of an existing matrix, this particular example is probably not the one to look at.

Hi mecej4,

Tianxiong_Lu — Fri, 18 Sep 2015 02:21:07 GMT

Hi mecej4,

Thank you for this information.

I have noticed the existence of L in comment, but missing the next line:

! L is the lower triangular part of matrix A with unit diagonal

and I can not find anything about matrix L in example code, so I misunderstand the true idea of example.

Now, it is work fine. Thanks again.

Best regards,

Tianxiong

Furthmore, I want solve the

Tianxiong_Lu — Fri, 18 Sep 2015 09:33:28 GMT

Furthmore, I want solve the full matrix equations rather than lower triangular matrix. And I just comment the matrix descriptor code and set its type to SPARSE_MATRIX_TYPE_GENERAL. But, unfortunately, the result is a wrong answer. What's the matter?

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "mkl_spblas.h"

int main() {
    //*******************************************************************************
    //     Declaration and initialization of parameters for sparse representation of
    //     the matrix A in the compressed sparse row format:
    //*******************************************************************************
#define M 5
#define N 5
#define NNZ 13
    //*******************************************************************************
    //    Sparse representation of the matrix A
    //*******************************************************************************
    double csrVal[NNZ]    = { 1.0, -1.0,     -3.0,
                             -2.0,  5.0,
                                         4.0, 6.0, 4.0,
                             -4.0,       2.0, 7.0,
                                    8.0,          -5.0 };
    MKL_INT    csrColInd[NNZ] = { 0,      1,        3,
                              0,      1,
                                           2,   3,   4,
                              0,           2,   3,
                                      1,             4 };
    MKL_INT    csrRowPtr[M+1] = { 0, 3, 5, 8, 11, 13 };
    MKL_INT    csrRowEndPtr = { 2, 4, 7, 10, 12 };
    // Descriptor of main sparse matrix properties
    struct matrix_descr descrA;
    // Structure with sparse matrix stored in CSR format
    sparse_matrix_t       csrA;
    //*******************************************************************************
    //    Declaration of local variables:
    //*******************************************************************************
    double x  = { 1.0, 5.0, 1.0, 4.0, 1.0};
    double y  = { 0.0, 0.0, 0.0, 0.0, 0.0};
    double y1  = { 0.0, 0.0, 0.0, 0.0, 0.0};
    double alpha = 1.0, beta = 0.0, error = 0.0;
    MKL_INT    i;

    printf( "\n EXAMPLE PROGRAM FOR mkl_sparse_d_trsv \n" );
    printf( "---------------------------------------------------\n" );
    printf( "\n" );
    printf( "   INPUT DATA FOR mkl_sparse_d_trsv    \n" );
    printf( "   LOWER TRIANGULAR SPARSE MATRIX \n" );
    printf( "   WITH UNIT DIAGONAL             \n" );
    printf( "   ALPHA = %4.1f  BETA = %4.1f    \n", alpha, beta );
    printf( "   SPARSE_OPERATION_NON_TRANSPOSE \n" );
    printf( "   Input vector                   \n" );
    for ( i = 0; i < N; i++ )
    {
        printf( "%7.1f\n", x );
    };

    // Create handle with matrix stored in CSR format
    mkl_sparse_d_create_csr ( &csrA, SPARSE_INDEX_BASE_ZERO,
                                    N,  // number of rows
                                    M,  // number of cols
                                    csrRowPtr,
                                    csrRowPtr+1,
                                    csrColInd,
                                    csrVal );

    // Create matrix descriptor
    //descrA.type = SPARSE_MATRIX_TYPE_TRIANGULAR;
    //descrA.mode = SPARSE_FILL_MODE_LOWER;
    //descrA.diag = SPARSE_DIAG_UNIT;
    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;

    // Compute y = alpha * A^{-1} * x
    mkl_sparse_d_trsv ( SPARSE_OPERATION_NON_TRANSPOSE,
                        alpha,
                        csrA,
                        descrA,
                        x,
                        y );

    printf( "                                   \n" );
    printf( "   OUTPUT DATA FOR mkl_sparse_d_trsv \n" );

    // y should be equal { 1.0, 7.0, 1.0, 6.0, -55.0 }
    for ( i = 0; i < N; i++ )
    {
        printf( "%.8f\n", y );
    };

    // For validation perform y1 = A * y
    mkl_sparse_d_mv ( SPARSE_OPERATION_NON_TRANSPOSE,
                      alpha,
                      csrA,
                      descrA,
                      y,
                      beta,
                      y1 );

    printf( "                                   \n" );
    printf( "   OUTPUT DATA FOR mkl_sparse_d_mv \n" );

    // y should be equal { 1.0, 7.0, 1.0, 6.0, -55.0 }
    for ( i = 0; i < N; i++ )
    {
        printf( "%7.1f\n", y1 );
    };

    // Release matrix handle and deallocate matrix
    mkl_sparse_destroy ( csrA );

    for ( i = 0; i < N; i++ )
    {
        error += x-y1;
    };
    if (error > 1e-8)
    {
        printf("\n   VALIDATION FAILED\n");
    } else
    {
        printf("\n   VALIDATION PASSED\n");
    }
    printf( "---------------------------------------------------\n" );
    return 0;
}

It give me an unexpected result:

   Input vector
    1.0
    5.0
    1.0
    4.0
    1.0

   OUTPUT DATA FOR mkl_sparse_d_trsv
18.00000000
5.00000000
-27.00000000
4.00000000
1.00000000

   OUTPUT DATA FOR mkl_sparse_d_mv
    1.0
  -11.0
  -80.0
  -98.0
   35.0

   VALIDATION FAILED
-------------------------------------------

Correct answer should be

-2.00148809523810
0.199404761904762
1.73139880952381
-1.06696428571429
0.119047619047619

As far as I can tell, mkl

mecej4 — Fri, 18 Sep 2015 15:13:28 GMT

As far as I can tell, mkl_sparse_?_trsv() and mkl_sparse_?_trmm() are new in MKL 11.3. The documentation pages for some of these new functions, I feel, need improvement. The "tr" in the function names is significant, standing for "triangular", as is conventional in naming Lapack functions/subroutines.

My guess is that, in those cases where the operation "op" involves matrix inversion (conceptually, at least), the operand must be triangular (or, as a special case, diagonal). These functions cannot solve simultaneous sparse linear equations when the matrix is not triangular. For solving sparse linear equations in the general case, when the matrix is square and symmetric or non-symmetric, positive-definite or not, you have to use Pardiso (directly or through the DSS interface).

Intel, please confirm or refute, and please consider revising MKL manual pages (for example, https://software.intel.com/en-us/node/590041 ) to make this restriction clear. If my conjecture is wrong, we need to be given some guidance on when to use Pardiso and when to use mkl_sparse_?_trsv, and a statement as to why the name contains "tr", in contravention of the normal Lapack naming schemes.

Hi Tianxiong, Mecej4,

Ying_H_Intel — Mon, 21 Sep 2015 02:08:58 GMT

Hi Tianxiong, Mecej4,

You are right, the tr functions in SPMV 2, BLAS 2 -3 are only valid for triangular system.

sv solving a single triangular system (Level 2)

sm solving triangular systems with multiple right-hand sides (Level 3).

So you can't use it for full sparse matrix. I have asked mkl docs team to add such discription under such functions. Thanks a lot for the raising the question,.

If you want to solve the full matrix equations, you may need to call lapack routines for dense matrix or sparse solver like PARDISO, DSS etc.

Best Regards,
Ying