#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'*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 <n ; i++){
        for(j = 0; j <m ; j++)
         printf("C[%d][%d]=%f  ", i,j,A[i*n+j]);
        printf("\n");
    }
}