Индексация на базе Intel MKL SpareBlas mm CSR не работает

Я тестирую функции Intel MKL в тестовой программе C и обнаружил, что просто не могу сделать spareblas: mkl_scsrmm функцией CSR работа индексирования на основе одной. Я использую CSR с вариантами val, columns, pntrb и pntre. Оригинальные примеры были помещены в:

...mkl\examples\examples_core_c\spblasc\source\cspblas_scsr.c

Это первый код для индексации с отсчетом от нуля:

Пример №1

#include <stdio.h>
#include "mkl_types.h"
#include "mkl_spblas.h"

int main() {
    #define M 2        
    #define NNZ 4
        
        MKL_INT        m = M, nnz = NNZ;
        
        float        values[NNZ]     = {2.0,4.0,4.0,2.0};
        MKL_INT        columns[NNZ]  = {1,2,1,2};
        MKL_INT        rowIndex[M+1] = {1,3,5};
        
    #define N 2      
      
        MKL_INT        n = N;
        float         b[M][N]    = {2.0, 1.0, 5.0, 2.0};
        float         c[M][N]    = {0.0, 0.0, 0.0, 0.0};
        float         alpha = 1.0, beta = 0.0;
        char          transa, uplo, nonunit;
        char          matdescra[6];
        MKL_INT       i, j, is;
        
    transa = 'N';
    
        matdescra[0] = 'S';
        matdescra[1] = 'L';
        matdescra[2] = 'N';
        matdescra[3] = 'F';
        
        mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, values, columns, rowIndex, &(rowIndex[1]), &(b[0][0]), &n,  &beta, &(c[0][0]), &n);

        printf("                             \n");
        printf("   OUTPUT DATA FOR MKL_SCSRMM\n");
        for (i = 0; i < m; i++) {
            for (j = 0; j < n; j++) {
                printf("%7.1f", c[i][j]);
            };
            printf("\n");
        };
        return 0;
}

Результаты, которые я получаю, таковы:

Индексирование с отсчетом от нуля (правильный вариант):

24.0 10.0

18.0 8.0

Единая индексация:

8.0 10.0

18.0 24.0

Хотя кажется, что это только меняет положение диагональных элементов, с матрицей 3x3 решение полностью отличается от правильного. Я подозревал, что это может быть что-то с форматом ввода матрицы b. Я думаю, что не хватает ясности в описании массива b для функции mkl_scsrmm, размещенном в справочнике по MKL. Таким образом, я изменил формат b в этом примере, и это сработало (я расположил элементы в следующем порядке: `{2.0, 5.0, 1.0, 2.0}). Но я сделал то же самое для другого примера 3x3, который я закодировал, и это не сработало. Я думаю, что это может быть просто совпадением. Я действительно не знаю, что делать с этой проблемой, я хотел бы понять, что здесь происходит.

Рекомендации:

Формат CSR

https://software.intel.com/en-us/node/471374

Запасная функция Blas mkl_scsrmm

https://software.intel.com/sites/products/documentation/doclib/iss/2013/mkl/mklman/hh_goto.htm#GUID-78C55D9B-86FF-4A9F-B5D5-D2F61B9314FC.htm< /а>

Запасные аспекты интерфейса Blas

Отн. ="nofollow noreferrer">https://software.intel.com/sites/products/documentation/doclib/iss/2013/mkl/mklman/hh_goto.htm#GUID-34C8DB79-0139-46E0-8B53-99F3BEE7B2D4.htm< /а>

А вот еще один пример, 3х3:

Пример #2

//  matrix A
//
//  2   4   3  
//  4   2   1
//  3   1   6
//
//  matrix B
//
//  2   1   3
//  4   5   6
//  7   8   9
//  
//  ZERO-BASED INDEXING
//  
//  a = {2 4 3 4 2 1 3 1 6}
//  columns= {0 1 2 0 1 2 0 1 2}
//  idexRow = {0 3 6 9}
//
//  b = {2 1 3 4 5 6 7 8 9} (row order array)
//
//  We print the array in row-major order
//
//  ONE-BASED INDEXING
//  
//  a = {2 4 3 4 2 1 3 1 6}
//  columns={1 2 3 1 2 3 1 2 3}
//  indexRow = {0 3 6 9}
//
//  b = {2 4 7 1 5 8 3 6 9} (column order array)
//  
//  We print the array in column-major order (because the resoult is in column major order, ie transposed)
//
//
//

#include <stdio.h>
#include "mkl_types.h"
#include "mkl_spblas.h"

int main() 
{

#define M 3        
#define NNZ 9 
#define N 3 

        MKL_INT     m = M, nnz = NNZ, n=N;
        float       a[NNZ]    = {2.0,4.0,3.0,4.0,2.0,1.0,3.0,1.0,6.0};
        MKL_INT     columns[NNZ]  = {0,1,2,0,1,2,0,1,2};
        MKL_INT     rowIndex[M+1] = {0,3,6,9};

        float       b[M][N] = {2.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
        float       c[M][N] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

        float       alpha = 1.0, beta = 0.0;
        MKL_INT     i, j;
        char        transa;
        char        matdescra[6];

        float       a1[NNZ]   = {2.0,4.0,3.0,4.0,2.0,1.0,3.0,1.0,6.0};
        MKL_INT     columns1[NNZ]  = {1,2,3,1,2,3,1,2,3};
        MKL_INT     rowIndex1[M+1] = {1,4,7,10};
        float       b1[M][N]    = {2.0, 4.0, 7.0, 1.0, 5.0, 8.0, 3.0, 6.0, 9.0};
        float       c1[M][N]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

//********************************
        //ZERO-BASED INDEXING
//********************************
        
        transa = 'n';
        matdescra[0] = 's';
        matdescra[1] = 'l';
        matdescra[2] = 'n';
        matdescra[3] = 'c';

        mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, a, columns, rowIndex, &(rowIndex[1]), &(b[0][0]), &n,  &beta, &(c[0][0]), &n);

        printf("                             \n");
        printf("   Right Solution: ZERO-BASED: C    \n");
        for (i = 0; i < m; i++)
        {
            for (j = 0; j < n; j++) {
                printf("%7.1f", c[i][j]);
            };
            printf("\n");
        };

        printf("                             \n");
        printf("   ZERO-BASED: C'    \n");
        for (i = 0; i < m; i++) 
        {
            for (j = 0; j < n; j++) 
            {
                printf("%7.1f", c[j][i]);
            };
            printf("\n");
        };

//********************************
        //ONE-BASED INDEXING
//********************************
        
        matdescra[3] = 'f';

        mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, a1, columns1, rowIndex1, &(rowIndex1[1]), &(b1[0][0]), &n,  &beta, &(c1[0][0]), &n);

        printf("                             \n");
        printf("   ONE-BASED: C    \n");
        for (i = 0; i < m; i++) 
        {
            for (j = 0; j < n; j++) 
            {
                printf("%7.1f", c1[i][j]);
            };
            printf("\n");
        };

        printf("                             \n");
        printf("   ONE-BASED: C'    \n");
        for (i = 0; i < m; i++) 
        {
            for (j = 0; j < n; j++) 
            {
                printf("%7.1f", c1[j][i]);
            };
            printf("\n");
        };
        return 0;
}

person Adrian    schedule 23.10.2014    source источник


Ответы (1)


Я задал тот же вопрос на форуме Intel, и мне там помогли и я добрался до решения проблемы. Дело в том, что при вызове подпрограммы из интерфейса C с индексацией начиная с нуля вы можете отправить матрицу, хранящуюся в массиве, в построчном порядке (собственный C хранение массива), и когда вы вызываете подпрограмму с индексацией на основе единицы, вы должны хранить матрицу в порядке по столбцам. Это меняет способ хранения матриц B и C и способ сохранения результата. Для матрицы A меняется только индексация (с 0 на 1). Из документации Intel вы можете подумать, что интерфейс C всегда принимает порядок строк для обоих типов индексации.

Обратите внимание, что в целом упорядочение по столбцам не то же самое, что и хранение транспонированной матрицы в порядке по строкам< /strong> (то же самое, если матрицы квадратные).

person Adrian    schedule 26.10.2014
comment
Большое спасибо за публикацию, а также за ответ на свой вопрос. Вы мне очень помогли, я пару часов безуспешно пытался использовать хранилище по столбцам, читая документацию. Вы знаете, указано ли это где-нибудь в документации MKL? - person fcdimitr; 15.11.2017