《算法导论》——矩阵乘法的Strassen算法

前言：

　　很多朋友看到我写的《算法导论》系列，可能会觉得云里雾里，不知所云。这里我再次说明，本系列博文时配合《算法导论》一书，给出该书涉及的算法的c++实现。请结合《算法导论》一书阅读该系列博文。我这里有该书的电子版，有需要的朋友可以留言。

正题：

　　今天讨论的算法是矩阵乘法的Strassen算法，该算法的精髓在于减少n/2矩阵*n/2矩阵的次数。首先，作一些写该算法的基础工作：

/*
    * 矩阵的加法运算
    */
    void Add(int** matrixA, int** matrixB, int** matrixResult,int length)
    {
        for(int i = 0; i < length; i++) {
            for(int j = 0; j < length; j++) {
                matrixResult[i][j] = matrixA[i][j] + matrixB[i][j];
            }
        }
    }

    /*
    * 矩阵的减法运算
    */
    void Sub(int** matrixA, int** matrixB, int** matrixResult,int length)
    {
        for(int i = 0; i < length; i++) {
            for(int j = 0; j < length; j++) {
                matrixResult[i][j] = matrixA[i][j] - matrixB[i][j];
            }
        }
    }

    /*
    * 矩阵乘法
    */
    void Mul(int** matrixA, int** matrixB, int** matrixResult){
        for(int i = 0; i < 2; ++i) {
            for(int j = 0; j < 2; ++j) {
                matrixResult[i][j] = 0;
                for(int k = 0; k < 2; ++k) {
                    matrixResult[i][j] += matrixA[i][k] * matrixB[k][j];
                }
            }
        }
    }

接着进入核心部分：

void Strassen(int** matrixA, int** matrixB, int** matrixResult,int length)
    {
        int halfLength=length/2;
        int** a11=new int*[halfLength];
        int** a12=new int*[halfLength];
        int** a21=new int*[halfLength];
        int** a22=new int*[halfLength];

        int** b11=new int*[halfLength];
        int** b12=new int*[halfLength];
        int** b21=new int*[halfLength];
        int** b22=new int*[halfLength];

        int** s1=new int*[halfLength];
        int** s2=new int*[halfLength];
        int** s3=new int*[halfLength];
        int** s4=new int*[halfLength];
        int** s5=new int*[halfLength];
        int** s6=new int*[halfLength];
        int** s7=new int*[halfLength];

        int** matrixResult11=new int*[halfLength];
        int** matrixResult12=new int*[halfLength];
        int** matrixResult21=new int*[halfLength];
        int** matrixResult22=new int*[halfLength];        

        int** temp=new int*[halfLength];
        int** temp1=new int*[halfLength];
        if(halfLength==1){
            Mul(matrixA, matrixB, matrixResult);
        }else{
            //首先将矩阵A，B 分为4块
            for(int i = 0; i < halfLength; i++) {    
                a11[i]=new int[halfLength];
                a12[i]=new int[halfLength];
                a21[i]=new int[halfLength];
                a22[i]=new int[halfLength];

                b11[i]=new int[halfLength];
                b12[i]=new int[halfLength];
                b21[i]=new int[halfLength];
                b22[i]=new int[halfLength];

                s1[i]=new int[halfLength];
                s2[i]=new int[halfLength];
                s3[i]=new int[halfLength];
                s4[i]=new int[halfLength];
                s5[i]=new int[halfLength];
                s6[i]=new int[halfLength];
                s7[i]=new int[halfLength];

                matrixResult11[i]=new int[halfLength];
                matrixResult12[i]=new int[halfLength];
                matrixResult21[i]=new int[halfLength];
                matrixResult22[i]=new int[halfLength];

                temp[i]=new int[halfLength];
                temp1[i]=new int[halfLength];
                for(int j = 0; j < halfLength; j++) {
                    a11[i][j]=matrixA[i][j];
                    a12[i][j]=matrixA[i][j+halfLength];
                    a21[i][j]=matrixA[i+halfLength][j];
                    a22[i][j]=matrixA[i+halfLength][j+halfLength];
                    b11[i][j]=matrixB[i][j];
                    b12[i][j]=matrixB[i][j+halfLength];
                    b21[i][j]=matrixB[i+halfLength][j];
                    b22[i][j]=matrixB[i+halfLength][j+halfLength];
                }                
            }

            //计算s1
            Sub(b12, b22, temp,halfLength);
            Strassen(a11, temp, s1,halfLength);
            //计算s2
            Add(a11, a12, temp,halfLength);
            Strassen(temp, b22, s2,halfLength);
            //计算s3
            Add(a21, a22, temp,halfLength);
            Strassen(temp, b11, s3,halfLength);
            //计算s4
            Sub(b21, b11, temp,halfLength);
            Strassen(a22, temp, s4,halfLength);
            //计算s5
            Add(a11, a22, temp1,halfLength);
            Add(b11, b22, temp,halfLength);
            Strassen(temp1, temp, s5,halfLength);
            //计算s6
            Sub(a12, a22, temp1,halfLength);
            Add(b21, b22, temp,halfLength);
            Strassen(temp1, temp, s6,halfLength);
            //计算s7
            Sub(a11, a21, temp1,halfLength);
            Add(b11, b12, temp,halfLength);
            Strassen(temp1, temp, s7,halfLength);

            //计算matrixResult11
            Add(s5, s4, temp1,halfLength);
            Sub(temp1, s2, temp,halfLength);
            Add(temp, s6, matrixResult11,halfLength);
            //计算matrixResult12
            Add(s1, s2, matrixResult12,halfLength);
            //计算matrixResult21
            Add(s3, s4, matrixResult21,halfLength);
            //计算matrixResult22
            Add(s5, s1, temp1,halfLength);
            Sub(temp1, s3, temp,halfLength);
            Sub(temp, s7, matrixResult22,halfLength);

            //结果送回matrixResult中
            for(int i = 0; i < halfLength; i++) {
                for(int j = 0; j < halfLength; j++) {
                    matrixResult[i][j]=matrixResult11[i][j];
                    matrixResult[i][j+halfLength]=matrixResult12[i][j];
                    matrixResult[i+halfLength][j]=matrixResult21[i][j];
                    matrixResult[i+halfLength][j+halfLength]=matrixResult22[i][j];
                }
                delete(a11[i]);
                delete(a12[i]);
                delete(a21[i]);
                delete(a22[i]);

                delete(b11[i]);
                delete(b12[i]);
                delete(b21[i]);
                delete(b22[i]);

                delete(s1[i]);
                delete(s2[i]);
                delete(s3[i]);
                delete(s4[i]);
                delete(s5[i]);
                delete(s6[i]);
                delete(s7[i]);

                delete(matrixResult11[i]);
                delete(matrixResult12[i]);
                delete(matrixResult21[i]);
                delete(matrixResult22[i]);

                delete(temp[i]);
                delete(temp1[i]);
            }
            delete(a11);
            delete(a12);
            delete(a21);
            delete(a22);

            delete(b11);
            delete(b12);
            delete(b21);
            delete(b22);

            delete(s1);
            delete(s2);
            delete(s3);
            delete(s4);
            delete(s5);
            delete(s6);
            delete(s7);

            delete(matrixResult11);
            delete(matrixResult12);
            delete(matrixResult21);
            delete(matrixResult22);

            delete(temp);
            delete(temp1);
        }
    }

该算法看着或许有些冗长，几乎一半都在进行动态指针的初始化和删除。利用该算法计算矩阵乘的时间复杂度为θ（n^lg7）。

测试一下吧：

#include "stdafx.h"
#include <iostream>
#include "SquareMatrix.h"

using namespace std;
using namespace dksl;


//STRASSEN矩阵乘法算法


const int N=8; //常量N用来定义矩阵的大小
int _tmain(int argc, _TCHAR* argv[])
{
    int **a=new int*[4];
    int **b=new int*[4];
    int **c=new int*[4];
    for(int i=0;i<4;i++)
        {
            a[i]=new int[4];
            b[i]=new int[4];
            c[i]=new int[4];
        for(int j=0;j<4;j++)
        {
            a[i][j]=1;
            b[i][j]=2;
        }
    }
    Strassen(a,b,c,4);
    for(int i=0;i<4;i++)
    {
        for(int j=0;j<4;j++)
            cout<<c[i][j]<<" ";
        cout<<endl;
    }
    system("PAUSE");
    return 0;
}