递归矩阵乘法c语言,递归矩阵乘法(Recursive matrix multiplication)

我读算法导论由CLRS。 本书展示了简单的分而治之矩阵乘法伪代码:

n = A.rows

let c be a new n x n matrix

if n == 1

c11 = a11 * b11

else partition A, B, and C

C11 = SquareMatrixMultiplyRecursive(A11, B11)

+ SquareMatrixMultiplyRecursive(A12, B21)

//...

return C

其中,例如,A11是尺寸为N / 2 XN / 2的子矩阵。 作者还暗示,我应该使用,而不是创造新的矩阵来表示子矩阵指数计算,所以我这样做:

#include

#include

template

struct Matrix

{

Matrix(size_t r, size_t c)

{

Data.resize(c, std::vector(r, 0));

}

void SetSubMatrix(const int r, const int c, const int n, const Matrix& A, const Matrix& B)

{

for(int _c=c; _c

{

for(int _r=r; _r

{

Data[_c][_r] = A.Data[_c][_r] + B.Data[_c][_r];

}

}

}

static Matrix SquareMultiplyRecursive(Matrix& A, Matrix& B, int ar, int ac, int br, int bc, int n)

{

Matrix C(n, n);

if(n == 1)

{

C.Data[0][0] = A.Data[ac][ar] * B.Data[bc][br];

}

else

{

C.SetSubMatrix(0, 0, n / 2,

SquareMultiplyRecursive(A, B, ar, ac, br, bc, n / 2),

SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc, n / 2));

C.SetSubMatrix(0, n / 2, n / 2,

SquareMultiplyRecursive(A, B, ar, ac, br, bc + (n / 2), n / 2),

SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));

C.SetSubMatrix(n / 2, 0, n / 2,

SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc, n / 2),

SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc, n / 2));

C.SetSubMatrix(n / 2, n / 2, n / 2,

SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc + (n / 2), n / 2),

SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));

}

return C;

}

void Print()

{

for(int c=0; c

{

for(int r=0; r

{

std::cout << Data[c][r] << " ";

}

std::cout << "\n";

}

std::cout << "\n";

}

std::vector<:vector> > Data;

};

int main()

{

Matrix A(2, 2);

Matrix B(2, 2);

A.Data[0][0] = 2;

A.Data[0][1] = 1;

A.Data[1][0] = 1;

A.Data[1][1] = 2;

B.Data[0][0] = 2;

B.Data[0][1] = 1;

B.Data[1][0] = 1;

B.Data[1][1] = 2;

A.Print();

B.Print();

Matrix C(Matrix::SquareMultiplyRecursive(A, B, 0, 0, 0, 0, 2));

C.Print();

}

它给了我不正确的结果,寿我不知道我在做什么错?

Answer 1:

// Recursive naive matrix multiplication in C, not strassen.

// 2013-Feb-15 Fri 12:28 moshahmed/at/gmail

#include

#include

#include

#include

#define M 2

#define N (1<

typedef int mat[N][N]; // mat[2**M,2**M] for divide and conquer mult.

typedef struct { int ra, rb, ca, cb; } corners; // for tracking rows and columns.

// set A[a] = k

void set(mat A, corners a, int k){

int i,j;

for(i=a.ra;i

for(j=a.ca;j

A[i][j] = k;

}

// set A[a] = [random(l..h)].

void randk(mat A, corners a, int l, int h){

int i,j;

for(i=a.ra;i

for(j=a.ca;j

A[i][j] = l + rand()% (h-l);

}

// Print A[a]

void print(mat A, corners a, char *name) {

int i,j;

printf("%s = {\n",name);

for(i=a.ra;i

for(j=a.ca;j

printf("%4d, ", A[i][j]);

printf("\n");

}

printf("}\n");

}

// Return 1/4 of the matrix: top/bottom , left/right.

void find_corners(corners a, int i, int j, corners *b) {

int rm = a.ra + (a.rb - a.ra)/2 ;

int cm = a.ca + (a.cb - a.ca)/2 ;

*b = a;

if (i==0) b->rb = rm; // top rows

else b->ra = rm; // bot rows

if (j==0) b->cb = cm; // left cols

else b->ca = cm; // right cols

}

// Naive Multiply: A[a] * B[b] => C[c], recursively.

void mul(mat A, mat B, mat C, corners a, corners b, corners c) {

corners aii[2][2], bii[2][2], cii[2][2];

int i, j, m, n, p;

// Check: A[m n] * B[n p] = C[m p]

m = a.rb - a.ra; assert(m==(c.rb-c.ra));

n = a.cb - a.ca; assert(n==(b.rb-b.ra));

p = b.cb - b.ca; assert(p==(c.cb-c.ca));

assert(m>0);

if (n==1) {

C[c.ra][c.ca] += A[a.ra][a.ca] * B[b.ra][b.ca];

return;

}

// Create the smaller matrices:

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

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

find_corners(a, i, j, &aii[i][j]);

find_corners(b, i, j, &bii[i][j]);

find_corners(c, i, j, &cii[i][j]);

}

}

// Now do the 8 sub matrix multiplications.

// C00 = A00*B00 + A01*B10

// C01 = A00*B01 + A01*B11

// C10 = A10*B00 + A11*B10

// C11 = A10*B01 + A11*B11

mul( A, B, C, aii[0][0], bii[0][0], cii[0][0] );

mul( A, B, C, aii[0][1], bii[1][0], cii[0][0] );

mul( A, B, C, aii[0][0], bii[0][1], cii[0][1] );

mul( A, B, C, aii[0][1], bii[1][1], cii[0][1] );

mul( A, B, C, aii[1][0], bii[0][0], cii[1][0] );

mul( A, B, C, aii[1][1], bii[1][0], cii[1][0] );

mul( A, B, C, aii[1][0], bii[0][1], cii[1][1] );

mul( A, B, C, aii[1][1], bii[1][1], cii[1][1] );

}

int main() {

mat A, B, C;

corners ai = {0,N,0,N};

corners bi = {0,N,0,N};

corners ci = {0,N,0,N};

//set(A,ai,2);

//set(B,bi,2);

srand(time(0));

randk(A,ai, 0, 2);

randk(B,bi, 0, 2);

set(C,ci,0); // set to zero before mult.

print(A, ai, "A");

print(B, bi, "B");

mul(A,B,C, ai, bi, ci);

print(C, ci, "C");

return 0;

}

Answer 2:

我找到了解决办法... SetSubMatrix是完全地不正确的:

void SetSubMatrix(const int r, const int c, const int rn, const int cn, const Matrix& A, const Matrix& B)

{

for(int _c=c; _c

{

for(int _r=r; _r

{

Data[_c][_r] = A.Data[_c-c][_r-r] + B.Data[_c-c][_r-r];

}

}

}

static Matrix SquareMultiplyRecursive(Matrix& A, Matrix& B, int ar, int ac, int br, int bc, int n)

{

Matrix C(n, n);

if(n == 1)

{

C.Data[0][0] = A.Data[ac][ar] * B.Data[bc][br];

}

else

{

C.SetSubMatrix(0, 0, n / 2, n / 2,

SquareMultiplyRecursive(A, B, ar, ac, br, bc, n / 2),

SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc, n / 2));

C.SetSubMatrix(0, n / 2, n / 2, n,

SquareMultiplyRecursive(A, B, ar, ac, br, bc + (n / 2), n / 2),

SquareMultiplyRecursive(A, B, ar, ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));

C.SetSubMatrix(n / 2, 0, n, n / 2,

SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc, n / 2),

SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc, n / 2));

C.SetSubMatrix(n / 2, n / 2, n, n,

SquareMultiplyRecursive(A, B, ar + (n / 2), ac, br, bc + (n / 2), n / 2),

SquareMultiplyRecursive(A, B, ar + (n / 2), ac + (n / 2), br + (n / 2), bc + (n / 2), n / 2));

}

return C;

}

Answer 3:

这里是我的答案,取决于递归矩阵乘法。 仅对于N = 2 ^ M,其中M> = 2

template <:size_t size>

int matrix_mul_recursive(int N, int i, int j, const int (&A)[size][size], const int (&B)[size][size], int (&C)[size][size]) {

if (N == 1) {

return *(const_cast(&(A[0][0])) + i) * (*(const_cast(&(B[0][0])) + j));

}

else {

const int H = N / 2;

const int T = (size * H);

int r = i / size;

int c = 0;

if (j < size) {

c = j;

}

else {

c = j % size;

}

C[r][c] += matrix_mul_recursive(H, i, j, A, B, C) +

matrix_mul_recursive(H, i + H, T + j, A, B, C);

C[r][c + H] += matrix_mul_recursive(H, i, j + H, A, B, C) +

matrix_mul_recursive(H, i + H, T + j + H, A, B, C);

C[r + H][c] += matrix_mul_recursive(H, T + i, j, A, B, C) +

matrix_mul_recursive(H, T + i + H, T + j, A, B, C);

C[r + H][c + H] += matrix_mul_recursive(H, T + i, j + H, A, B, C) +

matrix_mul_recursive(H, T + i + H, T + j + H, A, B, C);

}

return 0;

}

文章来源: Recursive matrix multiplication

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值