矩阵乘法 分块 c语言,转置矩阵的分块并行乘法（C语言实现），计算矩阵C[rawn][rawn]=A[rawm][rawn]'*B[rawm][rawn]，子块大小为S*T，其算法实现原理参加本代码的...

#include #define rawm 4

#define rawn 4

#define rawp 4

//子矩阵的大小为S行*T列

#define S 2 //块矩阵的行

#define T 4 //块矩阵的列

typedef float data_type;//定义为32位浮点数

//将矩阵定义为全局常量和变量

//data_type rawA[rawm][rawn] = {4,2,3,3,2,1,4,2};

//data_type rawB[rawm][rawn] = {3,2,3,3,2,1,1,4};

data_type rawC[rawn][rawn]={{0.0}}; //乘法：C=A*B

//将矩阵定义为全局常量和变量

//data_type rawA[rawm][rawn] = {4,2,3,3,2,1,4,2};

//data_type rawB[rawm][rawn] = {3,2,3,3,2,1,1,4};

/*方阵*/

data_type rawA[rawm][rawn] = { 3, 2, 2, 3, 3, 1, 3, 3, 4, 4, 3, 3, 2, 3, 3, 4 };

data_type rawB[rawm][rawn] = { 1, 2, 3, 1, 1, 2, 4, 1, 1, 2, 4, 3, 1, 4, 4, 1 };

data_type rawCAOB[rawn][rawn]={{0.0}}; //乘法：C=A’*B

//函数声明，否则在函数定义之前调用会出现警告信息

void MatixPrin(data_type *A,int m,int n);

void MatixPros();

//分块后子矩阵的个数h*l,A矩阵分为S*T的子矩阵，B矩阵分为T*S的子矩阵

static int col_M = 1;

static int row_N = 1;

//矩阵A、B分块后，不能全分块时，最后一行和最后一列的子矩阵的大小

static int col_last = 1;

static int row_last = 1;

//矩阵分块处理：计算分块后子矩阵块的个数

void MatixPros(){

//将矩阵rawA[rawm][rawn]分为C_M*R_N个大小为S*T的子块，ceil(x)函数返回不小于x的最小整数

if (rawm % S == 0) {

col_M = rawm / S;

} else {

col_M = rawm / S + 1;

}

//AC_M = ceil((double) rawm / (double) S); //矩阵A分块后的行数

row_N = ceil((double) rawn / (double) T); //矩阵A分块后的列数，即矩阵B分块后的行数

col_last = rawm - (col_M - 1) * S;//最后一行

row_last = rawn - (row_N - 1) * T;//最后一列

}

//2.1 子矩阵乘法 C=A'*B

void SMblock_MultCAOB(int si,int sj,int sk,int subm,int subn,int subp) {

int i, j, k;

for (j = 0; j < subn; j++){ //列号

for (i = 0; i < subm; i++) { //行号

for (k = 0; k < subp; k++) { //并行

//printf("子块乘：C[%d][%d]+=A[%d][%d]*B[%d][%d] \n",sj * T + j,sk * S + k,si*S + i,sj * T + j,si * T + i,sk*S + k);

//C[j * p + k]+= A[i*m + j] * B[i * p + k]; //参考

rawCAOB[sj * T + j][sk * T + k] += rawA[si*S + i][sj * T + j]* rawB[si*S + i][sk * T + k];

}

}

}

}

//2.2 分块矩阵运算：调用乘法实现分块矩阵乘法 C=A'*B，区别在于分块调用时循环次数为N*N*M，而不是M*M*N

//@data_type A[M][N], data_type B[M][N], data_type C[N][N]

void Mult_blkCAOB(data_type *A, data_type *B, data_type *C) {

int i, j, k;

int count=0;//循环计数

//循环 的顺序可根据需要更换，不会影响计算的结果

for (j = 0; j < row_N; j++) {

for (i = 0; i < col_M; i++) {

for (k = 0; k < row_N; k++) {

printf("\t 第%d层循环： ",++count);

printf(" 分块乘法：C[%d][%d]+=A[%d][%d]*B[%d][%d] \n",j,k,i,j,i,k);

int mblk=S,nblk=T,pblk=T;//默认当前参与运算的两个子矩阵块的大小，必须每次循环重新赋初值

//计算当前子块的大小为mblk*nblk

if ((i == col_M - 1)) {

mblk = col_last;

}

if (j == row_N - 1) {

nblk = row_last;

}

if (k == row_N - 1) {

pblk = row_last;

}

//分块矩阵乘法C=A'*B

//SMblock_MultCAOB(i, j, k, mblk, nblk, mblk);

SMblock_MultCAOB(i, j, k, mblk, nblk,pblk);

}

}

}

printf("\t 输出C=A'*B的结果 \n");

MatixPrin(*rawCAOB,rawn,rawn);

}

int main(void) {

//暴力测试C[m][p]=A'[m][n]*B[p][n]

MatixPros();//矩阵分块处理

Mult_blkCAOB(*rawA, *rawB, *rawCAOB); //分块矩阵的乘法 C=A'*B

//暴力测试C[m][p]=A'[n][m]*B[n][p]，注意矩阵下标的对应关系发生了变化

//R_MultA(*rawA, *rawB, *rawC, rawn, rawm, rawn);

//printf("C=A'*B的暴力计算结果\n");

//MatixPrin(rawC[0],rawn,rawn);

return 0;

}

//转置矩阵的暴力乘法：C=A'*B,注意矩阵行号列号的变化

//@data_type A[n][m], data_type B[n][p], data_type C[m][p]

void R_MultA(data_type *A, data_type *B, data_type *C, int m, int n, int p) {

int i, j, k;

for (j = 0; j < m; j++) { //矩阵A转置前的列号，即转置后的行号

for (i = 0; i < n; i++) { //B的列号

for (k = 0; k < p; k++) { //并行

//printf("C[%d]=%f",j * p + k,C[j * p + k]);//检测部分数组的初始值不是0.0

//若不添加此句，部分C的初值不一定为0.

//if(i==0){

//C[j * p + k]=0.0;

//}

C[j * p + k]+= A[i*m + j] * B[i * p + k];

//printf("\t C[%d][%d]=A[%d][%d]*B[%d][%d]+=%f*%f=%f \n",j,k,i,j,i,k,A[i*m + j],B[i * p + k],C[j*p+k]);

}

}

}

}

// 计算结果打印输出函数

void MatixPrin(data_type *A,int m,int n){

int i,j;

for(i = 0; i