矩阵链乘最优化实现

矩阵链乘最优化算法(括号化算法)，关键要找到A[i]...A[j]矩阵链相乘做最少乘法次数（存在m[][]中）的相乘顺序，记录在矩阵s[][]中。再利用递归定义矩阵链乘算法。递归的出口是只有一个矩阵（直接返回）或者两个矩阵（返回相乘后的结果矩阵）的情况。

//#include"OptimalMatrixMultiplication.h"
#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<fstream.h>

//using namespace std;

#define N 6  //N Matrices
#define MAX 15 // Biggest number of elements in row/column of matrices
typedef long (*Mat)[MAX];       ////////////!!! define a 2dim array pointer
int num;
Mat pp
; // To store the result matrix of any two matrices' multiplication
long A[N+1][MAX][MAX]; // store matrices to be calculated, A[0] is not used !!!

//find m[][] and s[][]
void Matrix_Chain_Order(int *p,int m[N+1][N+1],int s[N+1][N+1])
{
for (int i=1;i<=N;i++)
{
m[i][i] = 0;
}
for(int L=2;L<=N;L++)
{
for(i=1;i<=N-L+1;i++)
{
int j = i+L-1;
m[i][j] = (1<<30);            ////m[i][j] can't be bigger than 2^31-1  !!!
for(int k=i;k<j;k++)        /// There is no " ^ " symbol in C/C++ language  !!!
{                            //////////// pow(x,y) -> x^y
int q = m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
if(q<m[i][j])
{
m[i][j] = q;
s[i][j] = k;
}
}
}
}

}

//print the multiplication order
void Print_OptimalParens(int s[N+1][N+1],int i,int j)
{
if(i == j)
{
cout<<"A"<<i;
}
else
{
cout<<"(";
Print_OptimalParens(s,i,s[i][j]);
Print_OptimalParens(s,s[i][j]+1,j);
cout<<")";
}
}

// To generate matrices, all matrices have "the same" size
void Init_Matrix(long A[][MAX][MAX],int *p)
{
//read metrix data from file "MetrixData" generated by rand();
ifstream fin("MetrixData.txt");
//vector<int> vec;
long idata;
fin>>idata;
for(int i=1;i<=N;i++)
{
for(int j=0;j<p[i-1];j++) //matrices have different size!
{
for(int k=0;k<p[i];k++)
{
A[i][j][k] = idata;
fin>>idata;
}
}
}
}

Mat Matrix_Multiplication(Mat a,Mat b,Mat c)
{
for(int i=0;i<MAX;i++)
{
for(int j=0;j<MAX;j++)
{
for(int k=0;k<MAX;k++)
{
c[i][j]+=a[i][k]*b[k][j];
}
}
}
return c;
}

Mat Matrix_Chain_Multiplication(int i,int j,int s[N+1][N+1])
{
if(i == j)
{
return A[i];
}
else if(i+1 == j)
{
pp[num] = (Mat)malloc(MAX*MAX*sizeof(long));
if(pp[num] == NULL)
{
cout<<"Can't allocate memory!";
exit(1);
}
memset(pp[num],0,MAX*MAX*sizeof(long));
return Matrix_Multiplication(A[i],A[j],pp[num ++]);
// num ++; // This line can't be excuted !!!
}
else
{
Mat l_Mat = Matrix_Chain_Multiplication(i,s[i][j],s);
Mat r_Mat = Matrix_Chain_Multiplication(s[i][j]+1,j,s);

pp[num] = (Mat)malloc(MAX*MAX*sizeof(long));
if(pp[num] == NULL)
{
cout<<"Can't allocate memory!";
exit(1);
}
memset(pp[num],0,MAX*MAX*sizeof(long));

return Matrix_Multiplication(l_Mat,r_Mat,pp[num ++]);
}
}

//print 2_Dim array
void Print_Array(long b[][MAX])
{
for(int i=0;i<MAX;i++)
{
for(int j=0;j<MAX;j++)
{
cout<<b[i][j]<<" ";
}
cout<<endl;
}
}

void main()
{
int m[N+1][N+1];
int s[N+1][N+1];
int p[N+1]={10,8,5,7,6,12,9}; //store the size of matrices
Matrix_Chain_Order(p,m,s);

// print m[][]
for(int i=1;i<=N;i++)
{
for(int j=i+1;j<=N;j++)
{
cout<<" m["<<i<<"]["<<j<<"]: "<<m[i][j];
}
cout<<endl;
}
//print s[][]
for(i=1;i<=N;i++)
{
for(int j=i+1;j<=N;j++)
{
cout<<" s["<<i<<"]["<<j<<"]: "<<s[i][j];
}
cout<<endl;
}

Print_OptimalParens(s,1,6);
cout<<endl;

Init_Matrix(A,p);
//Print_Array(A[1]);
Matrix_Chain_Multiplication(1,6,s);
Print_Array(pp[--num]);

for(i=0;i<=num;i++)
{
free(pp[i]);
}

}

OptimalMatrixMultiplication.cpp