转置矩阵的分块并行乘法(C语言实现),计算矩阵C[rawn][rawn]=A[rawm][rawn]'*B[rawm][rawn],子块大小为S*T,其算法实现原理参加本代码的附件。
2016-06-21 00:38
856 查看
#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, af71 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] , data_type B[M] , data_type C 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] *B[p][n ] MatixPros();//矩阵分块处理 Mult_blkCAOB(*rawA, *rawB, *rawCAOB); //分块矩阵的乘法 C=A'*B //暴力测试C[m][p]=A' [m]*B [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 [m], data_type B [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"); } }
相关文章推荐
- 矩阵基础知识
- 转置矩阵
- 023 A转置矩阵=A的性质(三大性质)
- 矩阵及其压缩存储
- 【Numpy学习记录】np.transpose讲解
- [leetcode] 【分治法】 50. Pow(x, n)
- C/C++代码命名和格式规范
- PAT乙级练习题B1013.数素数
- c++ String 转 char*
- c++ 深入理解虚函数
- 关于C++内联函数
- C语言的一些杂货
- 关于C++内联函数
- C++对C的增强
- C++Primer 第6章笔记整理
- c语言学习笔记32
- C++ 虚函数表解析
- C++中数组名、指针的引用传递
- 在Eclipse中运行C++程序出现"Launch failed. Binary not foud"和"Program file not Specified"的问题
- 批注:C++中复制构造函数与重载赋值操作符总结:默认浅拷贝,带指针的需要深拷贝