稀疏矩阵 part 3

▶ 各种稀疏矩阵数据结构下 y(n,1) = A(n,m) * x(m,1) 的实现，CPU版本

● MAT 乘法

int dotCPU(const MAT *a, const MAT *x, MAT *y)
{
    checkNULL(a); checkNULL(x); checkNULL(y);
    if (a->col != x->row)
    {
        printf("dotMATCPU dimension mismatch!
");
        return 1;
    }
    
    y->row = a->row;
    y->col = x->col;
    for (int i = 0; i < a->row; i++)
    {
        format sum = 0;
        for (int j = 0; j < a->col; j++)        
            sum += a->data[i * a->col + j] * x->data[j];                
        y->data[i] = sum;        
    }
    COUNT_MAT(y);
    return 0;
}

● CSR 乘法

int dotCPU(const CSR *a, const MAT *x, MAT *y)
{
    checkNULL(a); checkNULL(x); checkNULL(y);
    if (a->col != x->row)
    {
        printf("dotCSRCPU dimension mismatch!
");
        return 1;
    }
    
    y->row = a->row;
    y->col = x->col;
    for (int i = 0; i < a->row; i++)                            // i 遍历 ptr，j 遍历行内数据，A 中为 0 的元素不参加乘法
    {
        format sum = 0;
        for (int j = a->ptr[i]; j < a->ptr[i + 1]; j++)
            sum += a->data[j] * x->data[a->index[j]];
        y->data[i] = sum;
    }
    COUNT_MAT(y);
    return 0;
}

● ELL 乘法

int dotCPU(const ELL *a, const MAT *x, MAT *y)      // CPU ELL乘法
{
    checkNULL(a); checkNULL(x); checkNULL(y);
    if (a->colOrigin != x->row)
    {
        printf("dotELLCPU dimension mismatch!
");
        return 1;
    }

    y->row = a->col;
    y->col = x->col;
    for (int i = 0; i<a->col; i++)
    {
        format sum = 0;
        for (int j = 0; j < a->row; j++)
        {
            int temp = a->index[j * a->col + i];
            if (temp < 0)                                   // 跳过无效元素
                continue;
            sum += a->data[j * a->col + i] * x->data[temp];
        }
        y->data[i] = sum;
    }
    COUNT_MAT(y);
    return 0;
}

● COO 乘法

int dotCPU(const COO *a, const MAT *x, MAT *y)
{
    checkNULL(a); checkNULL(x); checkNULL(y);
    if (a->col != x->row)
    {
        printf("dotCOOCPU null!
");
        return 1;
    }

    y->row = a->row;
    y->col = x->col;
    for (int i = 0; i<a->count; i++)
        y->data[a->rowIndex[i]] += a->data[i] * x->data[a->colIndex[i]];
    COUNT_MAT(y);
    return 0;
}

● DIA 乘法

int dotCPU(const DIA *a, const MAT *x, MAT *y)
{
    checkNULL(a); checkNULL(x); checkNULL(y);
    if (a->colOrigin != x->row)
    {
        printf("dotDIACPU null!
");
        return 1;
    }    
    y->row = a->row;
    y->col = x->col;
    int * inverseIndex = (int *)malloc(sizeof(int) * a->col);
    for (int i = 0, j = 0; i < a->row + a->col - 1; i++)
    {
        if (a->index[i] == 1)
        {
            inverseIndex[j] = i;
            j++;
        }
    }
    for (int i = 0; i < a->row; i++)
    {
        format sum = 0;
        for (int j = 0; j < a->col; j++)
        {
            if (i < a->row - 1 - inverseIndex[j] || i > inverseIndex[a->col - 1] - inverseIndex[j])
                continue;
            sum += a->data[i * a->col + j] * x->data[i + inverseIndex[j] - a->row + 1];
        }
        y->data[i] = sum;
    }
    COUNT_MAT(y);
    free(inverseIndex);
    return 0;
}