二维数组用vector容器，实现Strassen算法，效率很低！help。。

[横云断岭]

matrix.h,

/*
本文件定义了CMatrix类，使用STL模板容器 vector 来存储矩阵，可以实现任意 n 的Strassen算法矩阵相乘
当 n 为奇数时，调用普通矩阵算法。
example:
CMatrix A(8),B(*),C(8);
A.Trun_Random();
B.Trun_Random();
C = A*B;
cout<<C;
*/


#ifndef matrix_h_
#define matrix_h_
#include<vector>
#include<iostream>
#include  <time.h> 
using namespace std;


class CMatrix{
	//运算符重载
	friend ostream& operator << (ostream &out,CMatrix& M);
	friend CMatrix& operator + (CMatrix& M1,CMatrix& M2);
	friend CMatrix& operator - (CMatrix& M1,CMatrix& M2);
	friend CMatrix& operator * (CMatrix& M1,CMatrix& M2);  //Strassen算法矩阵相乘
public:
	int N;
	vector< vector<int> > A;  //二维
public:
	CMatrix(int n); //构造，生成n阶二维数组，元素值都为1
	CMatrix(int n,int a);  //构造，生成n阶二维数组，元素值都为a;
	~CMatrix(){ };
	CMatrix(CMatrix& M2);

	CMatrix& operator = (CMatrix& M2); //重载“=”

	void Trun_Random();   //使数组的元素全部变成随机数，0~9
	void Union(CMatrix& A11,CMatrix& A12,CMatrix& A21,CMatrix& A22);//联合4个小二维数组
	CMatrix& normalMutil(CMatrix& M2);   //普通矩阵乘法
};

#endif

[横云断岭]

#include"matrix.h"
#include<iostream>
#include "iomanip"

using namespace std;  

int getRandom() //得到随机数
{
	return rand();
}
void CMatrix::Trun_Random() //使矩阵内元素变为0~9的元素
{
	for(int i = 0;i<N;i++)
	{
		for(int j = 0;j<N;j++)
			A[i][j] = getRandom()%10;
	}
}

CMatrix::CMatrix(int n)
{
	N = n;
	A =  vector<vector<int>>(N,vector<int>(N,1));

}

CMatrix::CMatrix(int n, int a)
{
	N = n;
	A = vector<vector<int>>(N,vector<int>(N,a));
}
CMatrix::CMatrix(CMatrix& M2)
{
	N = M2.N;
//	A = vector<vector<int>>(N,vector<int>(N,1));
	vector<vector<int>>::iterator iter;
	for(iter = M2.A.begin();iter!=M2.A.end();iter++)
		A.push_back(*iter);
}

//重载 + 
CMatrix& operator + (CMatrix& M1,CMatrix& M2)
{
	CMatrix* temp = new CMatrix(M1.N);
	
	if(M1.N != M2.N)
		cout<<"N != M2.N"<<endl;

	temp->A = M2.A;
	
	for(int i = 0;i<M1.N;i++)
	{
		for(int j = 0;j<M1.N;j++)
			temp->A[i][j] += M1.A[i][j];
	}

	return *temp;
}
//重载 —
CMatrix& operator - (CMatrix& M1,CMatrix& M2)
{
	if(M1.N != M2.N)
		cout<<"M1.N != M2.N"<<endl;

	CMatrix* temp = new CMatrix(M1.N);
	
	temp->A = M1.A;
	
	for(int i = 0;i<M1.N;i++)
	{
		for(int j = 0;j<M1.N;j++)
			temp->A[i][j] -= M2.A[i][j];
	}

	return *temp;
}

//重载=
CMatrix& CMatrix::operator =(CMatrix &M2)
{
	if(N != M2.N)
		cout<<"N != M2.N"<<endl;
	A = M2.A;

	return *this;
}

//重载*
CMatrix& operator * (CMatrix& A,CMatrix& B)
{
	int N = A.N;
	if(A.N != B.N)
		cout<<"error!A.N != B.N"<<endl;

	if(N == 2) //当阶数为2时，做普通乘法运算
	{
		CMatrix* temp = new CMatrix(N);
		(*temp) = A.normalMutil(B);
		return *temp;
	}
	if(N%2 != 0) //当阶数为奇数时，做普通乘法运算
	{
		CMatrix* temp = new CMatrix(N);
		(*temp) = A.normalMutil(B);
		return *temp;
	}
	//M1~M7是中间结果
	CMatrix *M1 = new CMatrix(N/2),*M2 = new CMatrix(N/2),*M3 = new CMatrix(N/2),*M4 = new CMatrix(N/2),
		*M5 = new CMatrix(N/2),*M6 = new CMatrix(N/2),*M7 = new CMatrix(N/2);
	//将矩阵分为4块
	//A11,A12,A21,A22 是第一个矩阵 A 的分块
	//B11,B12,B21,B22 是第二个矩阵 B 的分块
	//C11,C12,C21,C22 是结果矩阵 C 的分块
	CMatrix *A11 = new CMatrix(N/2),*A12 = new CMatrix(N/2),*A21 = new CMatrix(N/2),*A22 = new CMatrix(N/2);
	CMatrix *B11 = new CMatrix(N/2),*B12 = new CMatrix(N/2),*B21 = new CMatrix(N/2),*B22 = new CMatrix(N/2);
	CMatrix *C11 = new CMatrix(N/2),*C12 = new CMatrix(N/2),*C21 = new CMatrix(N/2),*C22 = new CMatrix(N/2);
	CMatrix *C = new CMatrix(N);  //结果存放在C中

	//执行矩阵分块
	for(int i = 0;i<N/2;i++)
	{
		for(int j = 0;j<N/2;j++)
		{
			A11->A[i][j]=A.A[i][j];
			A12->A[i][j]=A.A[i][j+N/2];
			A21->A[i][j]=A.A[i+N/2][j];
			A22->A[i][j]=A.A[i+N/2][j+N/2];
			B11->A[i][j]=B.A[i][j];
			B12->A[i][j]=B.A[i][j+N/2];
			B21->A[i][j]=B.A[i+N/2][j];
			B22->A[i][j]=B.A[i+N/2][j+N/2];
		}
	}
	//计算M1~M7
	(*M1) = (*A11)*((*B12) - (*B22));
	(*M2) = ((*A11) +(*A12))*(*B22);
	(*M3 )= ((*A21) + (*A22))*(*B11);
	(*M4) = (*A22)*((*B21) - (*B11));
	(*M5) = ((*A11)+(*A22))*((*B11)+(*B22));
	(*M6) = ((*A12)-(*A22))*((*B21)+(*B22)); 
	(*M7) = ((*A11)-(*A21))*((*B11)+(*B12));

	(*C11) = (*M5)+(*M4)-(*M2)+(*M6);
	(*C12) = (*M1)+(*M2);
	(*C21) = (*M3)+(*M4);
	(*C22) = (*M5)+(*M1)-(*M3)-(*M7);

	//把C11,C12,C21,C22合并成结果矩阵C
	(*C).Union((*C11),(*C12),(*C21),(*C22));

	return *C;

}

void CMatrix::Union(CMatrix &A11, CMatrix &A12, CMatrix &A21, CMatrix &A22)
{//把四个小矩阵合成一个矩阵
	if(A11.N*2 != N)
		cout<<"error! A11.N*2 != N";
	for(int i = 0;i<N/2;i++)
	{
		for(int j = 0;j<N/2;j++)
		{
			A[i][j] = A11.A[i][j];
			A[i][j+N/2] = A12.A[i][j];
			A[i+N/2][j] = A21.A[i][j];
			A[i+N/2][j+N/2] = A22.A[i][j];
		}
	}
}
CMatrix& CMatrix::normalMutil(CMatrix &M2)  //普通矩阵乘法运算
{
	if(N != M2.N)
		cout<<"(N != M2.N!";

	CMatrix *re = new CMatrix(N,0);

	for(int i = 0;i<N;i++)
	{
		for(int j = 0;j<N;j++)
		{
			(*re).A[i][j] = 0;
			for(int k = 0;k<N;k++)
				(*re).A[i][j] += A[i][k]*M2.A[k][j];
		}
	}

	return *re;
}
ostream& operator <<(ostream& out,CMatrix& M)
{//重载〈〈，输出矩阵
	for(int i = 0;i<M.N;i++)
	{
		for(int j = 0;j<M.N;j++)
		{
		//	out.precision(4);
			out<<setw(4)<<M.A[i][j]<<" ";
		}
		cout<<endl;
	}
	return out;
}

[横云断岭]

上面的是matrix.cpp文件，下面的是main.cpp文件

/*
编译环境 studio 2008
最后修改 2009/11/4 
*/
#include"matrix.h"
#include<time.h>
#include<stdio.h>

int main()
{
	int n = 1;
	srand( (unsigned)time(NULL )); //设置随机数发生种子

	//clock_t start, end;

	//start = clock();

	//end = clock();
 //   printf("Pi = %lf\n用时为%.2f秒\n", pi, (float)(end-start)/CLOCKS_PER_SEC);


clock_t start, end;
start = clock();

	int a[100][100],b[100][100],c[100][100];

		for(int i = 0;i<100;i++)
	{
		for(int j = 0;j<100;j++)
		{
			a[i][j] = 0;
			for(int k = 0;k<100;k++)
				a[i][j] += b[i][k]*c[k][j];
		}
	}
end = clock();

printf("用时为%.2f秒\n", (float)(end-start)/CLOCKS_PER_SEC);



	cout<<"程序将演示用Strassen算法计算矩阵相乘。"<<endl;
	cout<<"程序可以计算任意阶的矩阵乘法，当n<=0，退出,请输入n：";
	cin>>n;

	while(n>0)
	{
	
		cout<<"当n = "<<n<<"时"<<endl;

		CMatrix A(n),B(n),C(n); //创建三个n*n的矩阵对象

		//A和B中的元素变成0~9的随机整数
		A.Trun_Random();
		B.Trun_Random();

		cout<<"随机生成的第一个矩阵是："<<endl;
	//	cout<<A;
		cout<<"随机生成的第二个矩阵是："<<endl;
	//	cout<<B;


clock_t start1, end1,start2,end2;

start1 = clock();

		C = A*B;   //Strassen算法计算矩阵相乘
end1 = clock();
		cout<<"Strassen算法计算结果："<<endl;

	//	cout<<C;

start2 = clock();
		C = A.normalMutil(B);  //普通乘法计算
end2 = clock();
		cout<<"普通乘法计算结果："<<endl;
	//	cout<<C;


printf("用时为%.6f秒\n",(float)(end1-start1)/CLOCKS_PER_SEC);

printf("用时为%.6f秒\n",(float)(end2-start2)/CLOCKS_PER_SEC);

		cout<<"程序可以计算任意阶的矩阵乘法，请输入n：";
		cin>>n;
	}
}

[jackyjkchen]

什么编译器，谈STL的效率别用VC6……

[横云断岭]

悲剧了，在main里实现了一个普通的100*100的矩阵乘法，结果0.1秒。。
下面是我用vector容器的一些测试结果：
在递归深度不高的情况下Strassen算法和传统算法时间比较：
N/T(s) 100 200 300 400 500 600 1000
Strassen算法 0.671 4.930 12.948 35.146 56.799 91.276 404.836
传统算法 0.577 4.555 15.538 37.238 72.884 125.986 610.492
递归深度 2 3 2 4 2 3 3
当递归深度很高时，Strassen算法表现很差：
N/T(s) 32 64 128
Strassen算法 0.811 5.772 40.684
传统算法 0.031 0.140 1.186
递归深度 5 6 7

[横云断岭]

[Quote=引用 3 楼 jackyjkchen 的回复:]
什么编译器，谈STL的效率别用VC6……
[/Quote]
studio 2008
CPU AMD5000+

[jackyjkchen]

[Quote=引用 5 楼 hengyunabc 的回复:]
引用 3 楼 jackyjkchen 的回复:
什么编译器，谈STL的效率别用VC6……

studio 2008
CPU AMD5000+
[/Quote]
这就不好说了，我过去也经常质疑标准库里的东西，老被人骂……

[ahao]

vector的随机存取效率应该还是很高的，
但你要注意的是，如果频繁改变2维vector的大小，那么效率会是一个噩梦
因为如果resize造成重新申请内存的话，会把里面所有的vector复制一遍

[横云断岭]

我自己仔细看了下，貌似是内存没有释放。。
再研究下。。

[yshuise]

Boost::uBLAs库
这个库的开发由一系列的类似的库来指导：

由Jack Dongarra等人开发的BLAS 库。
由Todd Veldhuizen开发的Blitz++ 库。
由Scott Haney 等人开发的POOMA 库。
由Jeremy Siek等人开发的MTL 库。

BLAS 似乎是基本的线性代数使用最广泛的库，所以它可以被称作de-facto standard。它的接口是面向过程的(procedural)，单独的函数是对一些基本的线性代数运算的抽象。由于它是使用Fortran语言实现以及它所进行的优化，BLAS似乎也是现在最快的库之一。因为我无决定以面向对象的方式来设计和实现我们的类库，所以技术方式上有显著的不同。然而，每个人都可以使用我们的库中的运算符表达所有的BLAS抽象并对效率进行比较。

Blitz++ 是一个使用C++实现的让人映象深刻的(impressive)库。它的主要设计似乎是面向多维数组以及包括张量(tensor)的相关的运算符。Blitz++库的作者声称由于它的实现技术使用表达式模板(expression template)和模板元编程技术(template metaprogram)，他的库的性能与对应的Fortan代码的性能相当或更好一些。然而我们有一些理由来开发一个我们自己设计和实现的方法。我们不知道是否有人使用Blitz++库来实现传统的线性代数运算。我们同样也假定由于Blitz++库所使用的实现风格(idioms)，即使在今天也需要最高级的C++编译器技术。另一方面，Blitz++库也使用我们相信，使用表达式模板技术(expression templates)是将抽象惩罚降低到一个可接受限度的必需的技术方法。

POOMA 的设计目标似乎是在许多部分对Blitz++库的大量部分进行并行计算。它从偏微分方程和理论物理领域中提取类来扩展Blitz++的概念。这种实现支持并行体系结构(parallel architectures)。

MTL 是另一个使用C++来支持基本的线性代数运算的类库。我们共同的观点是线性代数库应当提供与BLAS库相应的功能。另一方面，我们认为C++标准库的概念并没有支持所需要的数值计算的概念要求。另一个区别是MTL库当前并没有使用表达式模板技术(expression templates)。这可能导致两个问题中的一个：可能存在表现能力缺失或可能的性能损失。

[yshuise]

概览
原理
功能
矩阵和向量类型概览
矩阵和向量运算概览
有效的 uBLAS 和 进一步的信息
向量(Vector)
向量(Vector)
单位向量(Unit Vector)
零向量(Zero Vector)
标量向量(Scalar Vector)
稀疏向量(Sparse Vector)
映射向量(Mapped Vector)
压缩向量(Compressed Vector)
(Coordinate Vector)
向量代理(Vector Proxies)
向量范围(Vector Range)
向量切分(Vector Slice)
向量表达式(Vector Expressions)
向量表达式(Vector Expression)
向量引用(Vector References)
向量运算(Vector Operations)
向量简化(Vector Reductions)
矩阵(Matrix)
矩阵(Matrix)
方形矩阵(Identity Matrix)
零矩阵(Zero Matrix)
标量矩阵(Scalar Matrix)
三角矩阵(Triangular Matrix)
三角矩阵(Triangular Matrix)
三角矩阵适配器(Triangular Adaptor)
对称矩阵(Symmetric Matrix)
对称矩阵(Symmetric Matrix)
对称矩阵适配器(Symmetric Adaptor)
埃尔米特矩阵(Hermitian Matrix)
埃尔米特矩阵(Hermitian Matrix)
埃尔米特矩阵(Hermitian Adaptor)
带状矩阵(Banded Matrix)
带状矩阵(Banded Matrix)
带状矩阵(Banded Adaptor)
稀疏矩阵(Sparse Matrix)
映射矩阵(Mapped Matrix)
压缩矩阵(Compressed Matrix)
(Coordinate Matrix)
矩阵策略(Matrix Proxies)
矩阵行(Matrix Row)
矩阵列(Matrix Column)
向量范围(Vector Range)
向量切分(Vector Slice)
矩阵范围(Matrix Range)
矩阵切分(Matrix Slice)
矩阵表达式(Matrix Expressions)
矩阵表达式(Matrix Expression)
矩阵引用(Matrix References)
矩阵运算(Matrix Operations)
矩阵与向量运算(Matrix Vector Operations)
(矩阵与矩阵的运算Matrix Matrix Operations)
存储和特殊容器
无限数组(Unbounded Array)
有限数组(Bounded Array)
(范围Range)
切分(Slice)
稀疏存储(Sparse Storage)
缺省标准映射(Default Standard Map)
映射数组(Map Array)
运算和函数
特殊乘积(Special Products)
BLAS
uBLAS 概念定义
容器概念(Container Concepts)
向量(Vector)
矩阵(Matrix)
表达式概念(Expression Concepts)
标量表达式(Scalar Expression)
向量表达式(Vector Expression)
矩阵表达式(Matrix Expression)
存储概念(Storage Concept)
迭代器概念(Iterator Concepts)
(索引双向迭代器Indexed Bidirectional Iterator)
索引随机访问迭代器(Indexed Random Access Iterator)
索引双向列/行迭代器(Indexed Bidirectional Column/Row Iterator)
索引随机访问列/行迭代器(Indexed Random Access Column/Row Iterator)

[横云断岭]

[Quote=引用 11 楼 hityrj 的回复:]

您好……我也遇到这个问题了……请问你最后是怎么解决的……可以告诉我么……非常感谢……我用二维的vector进行普通的矩阵相乘……500*500的需要3分钟……慢死了
[/Quote]
话说忘记了。。
不过你可以考虑自已写一个操作二维数组的类。
可以考虑用boost::multi_array。

[hityrj]

您好……我也遇到这个问题了……请问你最后是怎么解决的……可以告诉我么……非常感谢……我用二维的vector进行普通的矩阵相乘……500*500的需要3分钟……慢死了