Strassen方法

最新推荐文章于 2024-03-11 20:36:50 发布

BigHome

最新推荐文章于 2024-03-11 20:36:50 发布

阅读量2.5k

点赞数

分类专栏：算法基础（C++）

本文链接：https://blog.csdn.net/qq_27016061/article/details/79578654

版权

算法基础（C++）专栏收录该内容

7 篇文章 0 订阅

订阅专栏

本人还没有写出Strassen的代码实现（蓝瘦香菇），在这里先给大家一个链接，至少目前来说是写的比较清楚地一个代码点击打开链接，一直在想在分解矩阵的时候通过下标计算实现，而不去创建12个新的n/2*n/2的矩阵，这样就不需要在花费O(n*n)时间复制矩阵元素，但是在递归的参数的传输一直有问题，暂时先搁置一下，在这里卡的时间太久了，如果小伙伴们有好的解决方法，可以一起交流一下。

上边提到的问题，想到了解决方法，重新开始介绍一下Strassen方法。

问题描述：

如果接触过矩阵的话，对矩阵的乘法一定不会陌生，这里边就不对矩阵的乘法进行解释了。计算一个n×n矩阵A和矩阵B的乘法，暴力的解决方法可以是直接用三个for循环解决，但这样会花费时间，我们现在需要找一个运行时间短些的方法；

解决方案：

我们直接看Strassen方法，这里我们使用了分治策略，至于其详细的解释可以看一看这个链接的内容点击打开链接，下边的代码也参考了一部分链接中的，但我们的解决方案是参考《算法导论》第四章分治策略，和上述链接中的代码的解决方案还是有些不同。具体的不同是本文的代码中，分解矩阵的时候通过下标计算实现，而不去创建12个新的n/2×n/2的矩阵。好了看解题步骤：

将输出矩阵分解为n/2×n/2的矩阵，输入矩阵A和B不需要分解，使用一个结构体来存储需要计算的下标区域（如上图所示的：A11 , A12，A21，A22）；
创建10个n/2×n/2的矩阵S1,S2,...,S10，存储内容如下：
S1=B12-B22 S3=A21+A22 S5=B21-B11 S7=A12-A22      S9=A11-A21
S2=A11+A12 S4=B21-B11 S6=B11+B22 S8=B21+B22      S10=B11+B12
通过之前的得到的矩阵递归的计算7个n/2×n/2矩阵积P1,P2,...,P7；
最后我们通过Pi矩阵的不同组合进行加减运算，计算出结果矩阵C的子矩阵C11 , C12，C21，C22。花费时间。

代码

#include<iostream>
using namespace std;
struct domain{
	int low_x;
	int high_x;
	int low_y;
	int high_y;	
};
template<typename T>
class Strassen_class{
public:
      void ADD(T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T** MatrixResult, int MatrixSize );
      void SUB(T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T** MatrixResult, int MatrixSize );
      void MUL( T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize );//朴素算法实现
      void FillMatrix( T** MatrixA, T** MatrixB, int length);//A,B矩阵赋值
      void PrintMatrix(T **MatrixA,int MatrixSize);//打印矩阵
      void Strassen(int N,T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T **MatrixC);//Strassen算法实现
};
template<typename T>
void Strassen_class<T>::ADD(T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T** MatrixResult, int MatrixSize )
{
    for ( int i = 0; i < MatrixSize; i++)
    {
        for ( int j = 0; j < MatrixSize; j++)
        {
            MatrixResult[i][j] =  MatrixA[domain_A.low_x+i][domain_A.low_y+j] + MatrixB[domain_B.low_x+i][domain_B.low_y+j];
        }
    }
}

template<typename T>
void Strassen_class<T>::SUB(T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T** MatrixResult, int MatrixSize )
{
    for ( int i = 0; i < MatrixSize; i++)
    {
        for ( int j = 0; j < MatrixSize; j++)
        {
           MatrixResult[i][j] =  MatrixA[domain_A.low_x+i][domain_A.low_y+j] + MatrixB[domain_B.low_x+i][domain_B.low_y+j];
        }
    }
}

template<typename T>
void Strassen_class<T>::MUL( T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize )
{
    for (int i=0;i<MatrixSize ;i++)
    {
        for (int j=0;j<MatrixSize ;j++)
        {
            MatrixResult[i][j]=0;
            for (int k=0;k<MatrixSize ;k++)
            {
                MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j];
            }
        }
    }
}
/*
c++使用二维数组，申请动态内存方法
申请
int **A;
A = new int *[desired_array_row];
for ( int i = 0; i < desired_array_row; i++)
     A[i] = new int [desired_column_size];

释放
for ( int i = 0; i < your_array_row; i++)
    delete [] A[i];
delete[] A;

*/
template<typename T>
void Strassen_class<T>::Strassen(int N,T** MatrixA,domain domain_A, T** MatrixB,domain domain_B, T **MatrixC)
{
	int HalfSize = N/2;
    int newSize = N/2;
	domain domain_1={0,HalfSize,0,HalfSize};  //分割数组中的第一组的界限
	domain domain_2={HalfSize,N,0,HalfSize};  //分割数组中的第二组的界限
	domain domain_3={0,HalfSize,HalfSize,N};  //分割数组中的第三组的界限
	domain domain_4={HalfSize,N,HalfSize,N};  //分割数组中的第四组的界限
    if ( N <= 64 )    //分治门槛，小于这个值时不再进行递归计算，而是采用常规矩阵计算方法
    {
        MUL(MatrixA,MatrixB,MatrixC,N);
    }
    else
    {
        T** S1;
        T** S2;
        T** S3;
        T** S4;
        T** S5;
        T** S6;
        T** S7;
        T** S8;        
        T** S9;
        T** S10;
        
        T** C11;
        T** C12;
        T** C21;
        T** C22;
        
        T** P1;
        T** P2;
        T** P3;
        T** P4;
        T** P5;
        T** P6;
        T** P7;
        T** AResult;
        T** BResult;

        //making a 1 diminsional pointer based array.
		S1 = new T *[newSize];
		S2 = new T *[newSize];
		S3 = new T *[newSize];
		S4 = new T *[newSize];
		S5 = new T *[newSize];
		S6 = new T *[newSize];
		S7 = new T *[newSize];
		S8 = new T *[newSize];        
		S9 = new T *[newSize];
		S10 = new T *[newSize];
		
		P1 = new T *[newSize];
		P2 = new T *[newSize];
		P3 = new T *[newSize];
		P4 = new T *[newSize];
		P5 = new T *[newSize];
		P6 = new T *[newSize];
		P7 = new T *[newSize];
        
        C11 = new T *[newSize];
        C12 = new T *[newSize];
        C21 = new T *[newSize];
        C22 = new T *[newSize];
 
        AResult = new T *[newSize];
        BResult = new T *[newSize];

        int newLength = newSize;

        //making that 1 diminsional pointer based array , a 2D pointer based array
        for ( int i = 0; i < newSize; i++)
        {
           	S1[i] = new T[newLength];
			S2[i] = new T[newLength];
			S3[i] = new T[newLength];
			S4[i] = new T[newLength];
			S5[i] = new T[newLength];
			S6[i] = new T[newLength];
			S7[i] = new T[newLength];
			S8[i] = new T[newLength];        
			S9[i] = new T[newLength];
			S10[i] = new T[newLength];
			
			P1[i] = new T[newLength];
			P2[i] = new T[newLength];
			P3[i] = new T[newLength];
			P4[i] = new T[newLength];
			P5[i] = new T[newLength];
			P6[i] = new T[newLength];
			P7[i] = new T[newLength];
	        
	        C11[i] = new T[newLength];
	        C12[i] = new T[newLength];
	        C21[i] = new T[newLength];
	        C22[i] = new T[newLength];
	 
	        AResult[i] = new T[newLength];
	        BResult[i] = new T[newLength];
        }
         //here we calculate P1..P7 matrices .
        //P1[][]
        SUB(MatrixB,domain_2,MatrixB,domain_4,S1,HalfSize);	//S1 
        Strassen(HalfSize,MatrixA,domain_1,S1,domain_1,P1); //P1=A11*S1

        //P2[][]
        ADD( MatrixA,domain_1,MatrixA,domain_2,S2, HalfSize); //S2
        Strassen(HalfSize,S2,domain_1,MatrixB,domain_4, P2);  //P2=S2*B22

        //P3[][]
        ADD( MatrixA,domain_3,MatrixA,domain_4,S3, HalfSize); //S3
		Strassen(HalfSize,S3,domain_1,MatrixB,domain_1, P3);  //P3=S3*B11
      
        //P4[][]
        SUB(MatrixB,domain_3,MatrixB,domain_1,S4,HalfSize);	//S4 
        Strassen(HalfSize,MatrixA,domain_4,S4,domain_1,P1); //P4=A22*S4
        
        //P5[][]
        ADD( MatrixA,domain_1,MatrixA,domain_4,S5, HalfSize); //S5
        ADD( MatrixB,domain_1,MatrixB,domain_4,S6, HalfSize); //S6
        Strassen(HalfSize,S5,domain_1,S6,domain_1,P5); //P5=S5*S6

        //P6[][]
        SUB( MatrixA,domain_2,MatrixA,domain_4,S7, HalfSize); //S7
        ADD( MatrixB,domain_3,MatrixB,domain_4,S8, HalfSize); //S8
        Strassen(HalfSize,S7,domain_1,S8,domain_1,P6); //P6=S7*S8
        
        //P7[][]
        SUB( MatrixA,domain_1,MatrixA,domain_3,S9, HalfSize); //S9
        ADD( MatrixB,domain_1,MatrixB,domain_2,S10, HalfSize); //S10
        Strassen(HalfSize,S9,domain_1,S10,domain_1,P6); //P7=S9*S10
        
        //C11 = P5 + P4 - (P2 - P6);
        ADD( P5,domain_1,P4,domain_1,AResult, HalfSize);
	    SUB( P2,domain_1,P6,domain_1,BResult, HalfSize);
        SUB( AResult,domain_1,BResult,domain_1,C11, HalfSize);
     

        //C12 = P1 + P2;
        ADD( P1,domain_1,P2,domain_1,C12,HalfSize);

        //C21 = P3 + P4;
        ADD( P3,domain_1,P4,domain_1,C21,HalfSize);

        //C22 = P5 + P1 - (P3 + P7);
        ADD( P5,domain_1,P1,domain_1,AResult,HalfSize);
        ADD( P3,domain_1,P7,domain_1,BResult,HalfSize);
		SUB( AResult,domain_1,BResult,domain_1,C22, HalfSize);
		
        //at this point , we have calculated the c11..c22 matrices, and now we are going to
        //put them together and make a unit matrix which would describe our resulting Matrix.
        //组合小矩阵到一个大矩阵
        for (int i = 0; i < N/2 ; i++)
        {
            for (int j = 0 ; j < N/2 ; j++)
            {
                MatrixC[i][j] = C11[i][j];
                MatrixC[i][j + N / 2] = C12[i][j];
                MatrixC[i + N / 2][j] = C21[i][j];
                MatrixC[i + N / 2][j + N / 2] = C22[i][j];
            }
        }

        // 释放矩阵内存空间
        for (int i = 0; i < newLength; i++)
        {
            delete[] S1[i];delete[] S2[i];delete[] S3[i];
            delete[] S4[i];delete[] S5[i];delete[] S6[i];
			delete[] S7[i];delete[] S8[i];delete[] S9[i];
			delete[] S10[i];
            delete[] C11[i];delete[] C12[i];delete[] C21[i];
            delete[] C22[i];
            delete[] P1[i];delete[] P2[i];delete[] P3[i];delete[] P4[i];
            delete[] P5[i];delete[] P6[i];delete[] P7[i];
            delete[] AResult[i];delete[] BResult[i] ;
        }
        delete[] S1;delete[] S2;delete[] S3;delete[] S4;
        delete[] S5;delete[] S6;delete[] S7;delete[] S8;
        delete[] S9;delete[] S10;
        delete[] C11;delete[] C12;delete[] C21;delete[] C22;
        delete[] P1;delete[] P2;delete[] P3;delete[] P4;delete[] P5;
        delete[] P6;delete[] P7;
        delete[] AResult;
        delete[] BResult ;

    }//end of else

}

template<typename T>
void Strassen_class<T>::FillMatrix( T** MatrixA, T** MatrixB, int length)
{
    for(int row = 0; row<length; row++)
    {
        for(int column = 0; column<length; column++)
        {

            MatrixB[row][column] = (MatrixA[row][column] = rand() %5);
            //matrix2[row][column] = rand() % 2;//ba hazfe in khat 50% afzayeshe soorat khahim dasht
        }

    }
}
template<typename T>
void Strassen_class<T>::PrintMatrix(T **MatrixA,int MatrixSize)
{
    cout<<endl;
    for(int row = 0; row<MatrixSize; row++)
    {
        for(int column = 0; column<MatrixSize; column++)
        {


            cout<<MatrixA[row][column]<<"\t";
            if ((column+1)%((MatrixSize)) == 0)
                cout<<endl;
        }

    }
    cout<<endl;
}
int main(int argc,char* argv[])
{
	Strassen_class<int> stra;//定义Strassen_class类对象
    int MatrixSize = 0;
	
    int** MatrixA;    //存放矩阵A
    int** MatrixB;    //存放矩阵B
    int** MatrixC;    //存放结果矩阵

    clock_t startTime_For_Normal_Multipilication ;
    clock_t endTime_For_Normal_Multipilication ;

    clock_t startTime_For_Strassen ;
    clock_t endTime_For_Strassen ;
    srand(time(0));

    cout<<"\n请输入矩阵大小(必须是2的幂指数值(例如:32,64,512,..): ";
    cin>>MatrixSize;

    int N = MatrixSize;//for readiblity.
	domain domain_main={0,N,0,N};
    //申请内存
    MatrixA = new int *[MatrixSize];
    MatrixB = new int *[MatrixSize];
    MatrixC = new int *[MatrixSize];

    for (int i = 0; i < MatrixSize; i++)
    {
        MatrixA[i] = new int [MatrixSize];
        MatrixB[i] = new int [MatrixSize];
        MatrixC[i] = new int [MatrixSize];
    }

    stra.FillMatrix(MatrixA,MatrixB,MatrixSize);  //矩阵赋值

  //*******************conventional multiplication test
        cout<<"朴素矩阵算法开始时钟:  "<< (startTime_For_Normal_Multipilication = clock());

        stra.MUL(MatrixA,MatrixB,MatrixC,MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)

        cout<<"\n朴素矩阵算法结束时钟: "<< (endTime_For_Normal_Multipilication = clock());

        cout<<"\n矩阵运算结果... \n";
        stra.PrintMatrix(MatrixC,MatrixSize);

  //*******************Strassen multiplication test
        cout<<"\nStrassen算法开始时钟: "<< (startTime_For_Strassen = clock());

        stra.Strassen( N, MatrixA, domain_main, MatrixB, domain_main, MatrixC ); //strassen矩阵相乘算法

        cout<<"\nStrassen算法结束时钟: "<<(endTime_For_Strassen = clock());


    cout<<"\n矩阵运算结果... \n";
    stra.PrintMatrix(MatrixC,MatrixSize);

    cout<<"矩阵大小 "<<MatrixSize;
    cout<<"\n朴素矩阵算法: "<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)<<" Clocks.."<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)/CLOCKS_PER_SEC<<" Sec";
    cout<<"\nStrassen算法:"<<(endTime_For_Strassen - startTime_For_Strassen)<<" Clocks.."<<(endTime_For_Strassen - startTime_For_Strassen)/CLOCKS_PER_SEC<<" Sec\n";
    system("Pause");

	return 0;
}

BigHome

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
Strassen方法

本人还没有写出Strassen的代码实现（蓝瘦香菇），在这里先给大家一个链接，至少目前来说是写的比较清楚地一个代码点击打开链接，一直在想在分解矩阵的时候通过下标计算实现，而不去创建12个新的n/2*n/2的矩阵，这样就不需要在花费O(n*n)时间复制矩阵元素，但是在递归的参数的传输一直有问题，暂时先搁置一下，在这里卡的时间太久了，如果小伙伴们有好的解决方法，可以一起交流一下。 ...
复制链接

扫一扫