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

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#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'*Bvoid 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*nblkif ((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 <n ; i++){        for(j = 0; j <m ; j++)         printf("C[%d][%d]=%f  ", i,j,A[i*n+j]);        printf("\n");    }}