稀疏矩阵的存储与计算 gaxpy

1, gaxpy 数学公式

 

其中：       ， ，

2, 具体实例

        

    

3，用稠密矩阵的方法

本节将用于验证第4节中的稀疏计算的结果

hello_gaxpy_dense.cpp

#include <stdio.h>
#include <stdlib.h>struct Matrix_SP
{float* val;   // 矩阵中非零元素，按列收集在一起；int* c;         // 某列非零元素开始的index in val[...]；int* r;          // val中每个元素的行号；int M;int N;
};void gen_sparse_matrix(int m, int n, float* A, int lda)
{//step1, randomfor(int i=0; i<m; i++){for(int j=0; j<n; j++){if((1.0*rand())/RAND_MAX < 0.6f)A[i + j*lda] = 0.0f;elseA[i + j*lda] = (1.0*rand())/RAND_MAX;}}
}void print_matrix(int m, int n, float* A, int lda)
{for(int i=0; i<m; i++){for(int j=0; j<n; j++){printf(" %7.4f", A[i + j*lda]);}printf("\n");}
}void gemm_cc(int M, int N, int K, float alpha, float*A, int lda, float* B, int ldb, float beta, float* C, int ldc)
{for(int i=0; i<M; i++){for(int j=0; j<N; j++){float sigma = 0.0f;for(int k=0; k<K; k++){sigma += A[i + k*lda]*B[k + j*ldb];}C[i + j*ldc] = alpha*sigma + beta*C[i + j*ldc];}}
}int main()
{int m = 6;int n = 5;int lda = m;// y and x are column vector;//y(6) = y(6) + A(6,5) * x(5)^t;float* y = NULL;float* x = NULL;float* A = NULL;A = (float*)malloc(lda*n*sizeof(float));y = (float*)malloc(m*sizeof(float));x = (float*)malloc(n*sizeof(float));//step 0, gen dense matrix;gen_sparse_matrix(m, n, A, lda);print_matrix(m, n, A, lda);//step 1, make dense matrix sparse;//step 2, gen y;//step 3, gen x;//step 4, dense gaxpy// 4.1 tmp = Ax;//4.2 y = y + tmp;//step 5, sparse gaxpy// 5.1 tmp = Ax// 5.2 y = y + tmp;free(A);free(y);free(x);return 0;
}

Makefile

all: g.outg.out: hello_dense_sparse_gaxpy.cppg++ -g $< -o $@.PHONY: clean
clean:-rm -rf g.out

4，用稀疏矩阵的方法

存储A

4.1 压缩列表示法

struct Matrix_SP
{float* val;   // 矩阵中非零元素，按列收集在一起；int* c;         // 某列非零元素开始的index in val[...]；int* r;          // val中每个元素的行号；int M;int N;
};

struct Matrix_SP 用于存储一个稀疏矩阵；

现在来用它存储 第2节 中的A矩阵；

4.2 稀疏 gemv 计算