码迷,mamicode.com
首页 > 编程语言 > 详细

【算法导论】矩阵乘法

时间:2015-02-11 12:45:16      阅读:265      评论:0      收藏:0      [点我收藏+]

标签:

离过年都不到十天了,还要等到这周五才能回家,想想也一年没回家了。从寒假开始到现在,已经有二十来天,这期间把2014年总结中的寒假计划也大多数完成了:The Element Of Style的阅读,三门数学课《随机过程》、《工程优化》、《数值分析》的算法实现。回家过年期间肯定不会写博客了,今天一看,这个月只写了三篇,于是乎今天必须再写一篇来完成这个月的基本工作量。言归正传,这篇文章写写选修课《算法设计》作业题中的矩阵乘法的三种方法。


矩阵乘法


  • 传统方法
    • 理论公式
      C=AB
      cij=k=1naikbkj
    • 算法实现
void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法,三重循环
{
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;//之所以每次调用都清零,是因为前面是循环调用,如果只调用一次就不需要
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}

  • 分块相乘法
    • 理论公式
      A=[A11A21A12A22],B=[B11B21B12B22]
      C=AB=[C11C21C12C22]
      C11=A11B11+A12B21C12=A11B12+A12B22C21=A21B11+A22B21C22=A21B12+A22B22
    • 算法实现
void BlockMatrix()//分块矩阵计算
{     
      for(int i=0;i<N/2;i++)              
         for(int j=0;j<N/2;j++)
            {
                A11[i][j]=A[i][j];
                A12[i][j]=A[i][j+N/2];
                A21[i][j]=A[i+N/2][j];
                A22[i][j]=A[i+N/2][j+N/2];
                B11[i][j]=B[i][j];
                B12[i][j]=B[i][j+N/2];
                B21[i][j]=B[i+N/2][j];
                B22[i][j]=B[i+N/2][j+N/2];

                C11[i][j]=0;
                C12[i][j]=0;
                C21[i][j]=0;
                C22[i][j]=0;
            }       //将矩阵A和B式分为四块

         MATRIX_Multiply(N/2,A11,B11, AA);
         MATRIX_Multiply(N/2,A12,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A11,B12, AA);
         MATRIX_Multiply(N/2,A12,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B11, AA);
         MATRIX_Multiply(N/2,A22,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B12, AA);
         MATRIX_Multiply(N/2,A22,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

    for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
        for(int j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]

    }

  • Strassen法
    • 理论公式
      M1=A11(B12?B22)M2=B22(A11+A12)M3=B11(A21+A22)M4=A22(B21?B11)M5=(A11+A22)(B11+B22)M6=(A12?A22)(B21+B22)M7=(A12?A21)(B11+B12)C11=M4+M5+M6?M2C12=M1+M2C21=M3+M4C22=M1+M5?M3?M7
    • 算法实现
void STRASSEN()  //STRASSEN函数
{
    int i,j;//,x;

    for(i=0;i<N/2;i++)              
        for(j=0;j<N/2;j++)
        {
            A11[i][j]=A[i][j];
            A12[i][j]=A[i][j+N/2];
            A21[i][j]=A[i+N/2][j];
            A22[i][j]=A[i+N/2][j+N/2];
            B11[i][j]=B[i][j];
            B12[i][j]=B[i][j+N/2];
            B21[i][j]=B[i+N/2][j];
            B22[i][j]=B[i+N/2][j+N/2];
        }       //将矩阵A和B式分为四块




    MATRIX_SUB(N/2,B12,B22,BB);         
    MATRIX_Multiply(N/2,A11,BB,M1);

    MATRIX_ADD(N/2,A11,A12,AA);
    MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

    MATRIX_ADD(N/2,A21,A22,AA);
    MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

    MATRIX_SUB(N/2,B21,B11,BB);
    MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

    MATRIX_ADD(N/2,A11,A22,AA);
    MATRIX_ADD(N/2,B11,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


    MATRIX_SUB(N/2,A12,A22,AA);
    MATRIX_ADD(N/2,B21,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

    MATRIX_SUB(N/2,A11,A21,AA);
    MATRIX_ADD(N/2,B11,B12,BB);
    MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
    //计算M1,M2,M3,M4,M5,M6,M7(递归部分)


    MATRIX_ADD(N/2,M5,M4,MM1);                
    MATRIX_SUB(N/2,M2,M6,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

    MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

    MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

    MATRIX_ADD(N/2,M5,M1,MM1);
    MATRIX_ADD(N/2,M3,M7,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

    for(i=0;i<N/2;i++)
        for(j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]



}

  • 完整程序实现
#include <iostream>
#include<ctime>

using namespace std;

const int N=32; //常量N用来定义方阵的大小
void output(int n,float C[][N]); //函数声明部分
void TraditionalMethod(float A[][N],float B[][N],float C[][N]);//传统的矩阵相乘
void BlockMatrix();
void STRASSEN();
void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2]);

float A[N][N];
float B[N][N];
float C[N][N];  //定义三个矩阵A,B,C

float A11[N/2][N/2],A12[N/2][N/2],A21[N/2][N/2],A22[N/2][N/2];
float B11[N/2][N/2],B12[N/2][N/2],B21[N/2][N/2],B22[N/2][N/2];
float C11[N/2][N/2],C12[N/2][N/2],C21[N/2][N/2],C22[N/2][N/2];
float M1[N/2][N/2],M2[N/2][N/2],M3[N/2][N/2],M4[N/2][N/2],M5[N/2][N/2],M6[N/2][N/2],M7[N/2][N/2];
float AA[N/2][N/2],BB[N/2][N/2],MM1[N/2][N/2],MM2[N/2][N/2];
void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]); //矩阵加法函数X+Y—>Z

void main()
{
  //初始化,使相乘的两个矩阵都为全1矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
        {
            A[i][j]=1;
            B[i][j]=1;
        }
    //将结果矩阵C初始化为全0矩阵
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;

    clock_t start, finish; //用于计时
    double   duration; 
    int loop=0;

    cout<<"当矩阵为"<<N<<"×"<<N<<",循环次数为10000时:"<<endl<<endl;

    //--------使用传统方法--------------//
    start = clock();  
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        TraditionalMethod(A,B,C);   //传统方法计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用传统方法:"<<endl;
    cout<<"所需时间为:"<<duration<<endl<<endl;

  //  output(N,C);  //输出计算结果

    //---------使用分块矩阵乘法------------//
    start = clock();  
    loop=0;
    while(loop<10000)//循环10000次,这里可以更改
    {
        loop++;
        BlockMatrix();   //分块矩阵计算

    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用分块相乘方法:"<<endl;
    cout<<"所需时间为:"<<duration<<endl<<endl;
    //  output(N,C);  //输出计算结果

    //-------使用strassen方法-------------------//
    start = clock();  
    loop=0;
    while(loop<10000)//当时间非常小时,需要加大循环次数,这里可以更改
    {
        loop++;
        STRASSEN();   //调用STRASSEN函数计算      
    }
    finish = clock(); 
    duration = (double)(finish - start) / CLOCKS_PER_SEC; 
    cout<<"使用strassen方法:"<<endl;
    cout<<"所需时间为:"<<duration<<endl;
//  output(N,C);  //输出计算结果
}

void TraditionalMethod(float A[][N],float B[][N],float C[][N])//传统方法,三重循环
{
    for(int i=0;i<N;i++)
        for(int j=0;j<N;j++)
            C[i][j]=0;//之所以每次调用都清零,是因为前面是循环调用,如果只调用一次就不需要
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}


void BlockMatrix()//分块矩阵计算
{     
      for(int i=0;i<N/2;i++)              
         for(int j=0;j<N/2;j++)
            {
                A11[i][j]=A[i][j];
                A12[i][j]=A[i][j+N/2];
                A21[i][j]=A[i+N/2][j];
                A22[i][j]=A[i+N/2][j+N/2];
                B11[i][j]=B[i][j];
                B12[i][j]=B[i][j+N/2];
                B21[i][j]=B[i+N/2][j];
                B22[i][j]=B[i+N/2][j+N/2];

                C11[i][j]=0;
                C12[i][j]=0;
                C21[i][j]=0;
                C22[i][j]=0;
            }       //将矩阵A和B式分为四块

         MATRIX_Multiply(N/2,A11,B11, AA);
         MATRIX_Multiply(N/2,A12,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C11); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A11,B12, AA);
         MATRIX_Multiply(N/2,A12,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C12); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B11, AA);
         MATRIX_Multiply(N/2,A22,B21, BB);
         MATRIX_ADD(N/2,AA,BB,C21); //矩阵加法函数X+Y—>Z

         MATRIX_Multiply(N/2,A21,B12, AA);
         MATRIX_Multiply(N/2,A22,B22, BB);
         MATRIX_ADD(N/2,AA,BB,C22); //矩阵加法函数X+Y—>Z

    for(int i=0;i<N/2;i++)//将上面计算得到的结果放入结果矩阵C中
        for(int j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]

    }



void output(int n,float C[][N]) //矩阵输出函数
{
    int i,j;
    cout<<"输出矩阵:"<<endl;
    for(i=0;i<n;i++)
    {
        cout<<endl;
        for(j=0;j<n;j++)
            cout<<C[i][j]<<"  ";
    }
    cout<<endl<<endl;

}


void MATRIX_Multiply(int n,float A[][N/2],float B[][N/2],float C[][N/2])//矩阵加法函数X*Y—>C
{


    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            C[i][j]=0;
            for(int k=0;k<n;k++)
            {
                C[i][j]=C[i][j]+A[i][k]*B[k][j];
            }
        }
    }

}



void MATRIX_ADD(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵加法函数X+Y—>Z
{
    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]+Y[i][j];
}

void MATRIX_SUB(int n,float X[][N/2],float Y[][N/2],float Z[][N/2]) //矩阵减法函数X-Y—>Z
{

    int i,j;
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            Z[i][j]=X[i][j]-Y[i][j];

}

void STRASSEN()  //STRASSEN函数
{


    int i,j;//,x;

    for(i=0;i<N/2;i++)              
        for(j=0;j<N/2;j++)
        {
            A11[i][j]=A[i][j];
            A12[i][j]=A[i][j+N/2];
            A21[i][j]=A[i+N/2][j];
            A22[i][j]=A[i+N/2][j+N/2];
            B11[i][j]=B[i][j];
            B12[i][j]=B[i][j+N/2];
            B21[i][j]=B[i+N/2][j];
            B22[i][j]=B[i+N/2][j+N/2];
        }       //将矩阵A和B式分为四块




    MATRIX_SUB(N/2,B12,B22,BB);         
    MATRIX_Multiply(N/2,A11,BB,M1);

    MATRIX_ADD(N/2,A11,A12,AA);
    MATRIX_Multiply(N/2,AA,B22,M2);//M2=(A11+A12)B22

    MATRIX_ADD(N/2,A21,A22,AA);
    MATRIX_Multiply(N/2,AA,B11,M3);//M3=(A21+A22)B11

    MATRIX_SUB(N/2,B21,B11,BB);
    MATRIX_Multiply(N/2,A22,BB,M4);//M4=A22(B21-B11)

    MATRIX_ADD(N/2,A11,A22,AA);
    MATRIX_ADD(N/2,B11,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M5);//M5=(A11+A22)(B11+B22)


    MATRIX_SUB(N/2,A12,A22,AA);
    MATRIX_ADD(N/2,B21,B22,BB);
    MATRIX_Multiply(N/2,AA,BB,M6);//M6=(A12-A22)(B21+B22)

    MATRIX_SUB(N/2,A11,A21,AA);
    MATRIX_ADD(N/2,B11,B12,BB);
    MATRIX_Multiply(N/2,AA,BB,M7);//M7=(A11-A21)(B11+B12)
    //计算M1,M2,M3,M4,M5,M6,M7(递归部分)


    MATRIX_ADD(N/2,M5,M4,MM1);                
    MATRIX_SUB(N/2,M2,M6,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C11);//C11=M5+M4-M2+M6

    MATRIX_ADD(N/2,M1,M2,C12);//C12=M1+M2

    MATRIX_ADD(N/2,M3,M4,C21);//C21=M3+M4

    MATRIX_ADD(N/2,M5,M1,MM1);
    MATRIX_ADD(N/2,M3,M7,MM2);
    MATRIX_SUB(N/2,MM1,MM2,C22);//C22=M5+M1-M3-M7

    for(i=0;i<N/2;i++)
        for(j=0;j<N/2;j++)
        {
            C[i][j]=C11[i][j];
            C[i][j+N/2]=C12[i][j];
            C[i+N/2][j]=C21[i][j];
            C[i+N/2][j+N/2]=C22[i][j];
        }                                            //计算结果送回C[N][N]   
}

运行结果如下图:
技术分享


原文:http://blog.csdn.net/tengweitw/article/details/43732253
作者:nineheadedbird

【算法导论】矩阵乘法

标签:

原文地址:http://blog.csdn.net/tengweitw/article/details/43732253

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!