算法导论C语言实现: 分治策略 -- 矩阵乘法的Strassen算法

先直接粘Code

4_2.c

#include <Windows.h>
#include <common.h>
#include "strassen.h"

void print_mutrix(const int *A, int r, int c)
{
	int i, j;
	printf("-----------------------------------\n");
	for (i=0; i < r; ++i) {
		for (j=0; j < c; ++j) {
			printf("%d\t", A[i*c + j]);
		}
		printf("\n");
	}
	printf("-----------------------------------\n");
}

//SQUARE-MATRIX-MULTIPLY(A,B)
void square_matrix_multiply(
	__in	int sm_size,
	__in	const int *sm_A,
	__in	const int *sm_B,
	__out	int *sm_C )
{
	int i, j, k;
	int p;
	for (i = 0; i < sm_size; ++i) {
		for (j = 0; j < sm_size; ++j) {
			//p = &sm_C[i*sm_size + j];
			p = 0;
			for (k = 0; k < sm_size; ++k) {
				p += (sm_A[i*sm_size + k]) * \
				      (sm_B[k*sm_size + j]);
			}
			sm_C[i*sm_size + j] += p;
		}
	}
}


// You must zero all bytes in sm_C before calling the recursive function.
//
//To avoid copying data, we define ROW_STEP:
//	Arow1			Arow2
// [----------------------][----------------------]
// .............
// 
//  A11row1     A12row1     A11row2     A12row2
// [----------][----------][----------][----------]
// |<--------row_step---->|
// .............

#define SM_SUB_AD11(ad, row_step, hs) \
	(ad)

#define SM_SUB_AD12(ad, row_step, hs) \
	((ad) + (hs))

#define SM_SUB_AD21(ad, row_step, hs) \
	((ad) + (row_step)*(hs))

#define SM_SUB_AD22(ad, row_step, hs) \
	((ad) + (row_step)*(hs) + (hs))

//SQUARE-MATRIX-MULTIPLY-RECURSIVE(A,B)
void square_matrix_multiply_recursive(
	__in	int sm_size,
	__in	int row_step,
	__in	const int *sm_A,
	__in	const int *sm_B,
	__out	int *sm_C )
{
	int hs = sm_size/2;
	if (sm_size == 1) {
		*sm_C += (*sm_A) * (*sm_B);
	} else {
		//C11 = F(A11, B11) + F(A12, B21)
		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD11(sm_A, row_step, hs),
			SM_SUB_AD11(sm_B, row_step, hs),
			SM_SUB_AD11(sm_C, row_step, hs));

		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD12(sm_A, row_step, hs),
			SM_SUB_AD21(sm_B, row_step, hs),
			SM_SUB_AD11(sm_C, row_step, hs));

		//C12 = F(A11, B12) + (F(A12, B22)
		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD11(sm_A, row_step, hs),
			SM_SUB_AD12(sm_B, row_step, hs),
			SM_SUB_AD12(sm_C, row_step, hs));

		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD12(sm_A, row_step, hs),
			SM_SUB_AD22(sm_B, row_step, hs),
			SM_SUB_AD12(sm_C, row_step, hs));

		//C21 = F(A21, B11) + F(A22, B21)
		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD21(sm_A, row_step, hs),
			SM_SUB_AD11(sm_B, row_step, hs),
			SM_SUB_AD21(sm_C, row_step, hs));

		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD22(sm_A, row_step, hs),
			SM_SUB_AD21(sm_B, row_step, hs),
			SM_SUB_AD21(sm_C, row_step, hs));

		//C22 = F(A21, B12) + F(A22, B22)
		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD21(sm_A, row_step, hs),
			SM_SUB_AD12(sm_B, row_step, hs),
			SM_SUB_AD22(sm_C, row_step, hs));

		square_matrix_multiply_recursive(
			hs, row_step,
			SM_SUB_AD22(