矩阵链乘最优化实现

Posted 2020-07-26 Y_SL

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

  1 //#include"OptimalMatrixMultiplication.h"
  2 #include<iostream.h>
  3 #include<stdio.h>
  4 #include<stdlib.h>
  5 #include<memory.h>
  6 #include<fstream.h>
  7 
  8 //using namespace std;
  9 
 10 #define N 6  //N Matrices
 11 #define MAX 15 // Biggest number of elements in row/column of matrices
 12 typedef long (*Mat)[MAX];       ////////////!!! define a 2dim array pointer
 13 int num;
 14 Mat pp[N]; // To store the result matrix of any two matrices\' multiplication
 15 long A[N+1][MAX][MAX]; // store matrices to be calculated, A[0] is not used !!!
 16 
 17 //find m[][] and s[][]
 18 void Matrix_Chain_Order(int *p,int m[N+1][N+1],int s[N+1][N+1])
 19 {
 20     for (int i=1;i<=N;i++)
 21     {
 22         m[i][i] = 0;
 23     }
 24     for(int L=2;L<=N;L++)
 25     {
 26         for(i=1;i<=N-L+1;i++)
 27         {
 28             int j = i+L-1;
 29             m[i][j] = (1<<30);            ////m[i][j] can\'t be bigger than 2^31-1  !!!
 30             for(int k=i;k<j;k++)        /// There is no " ^ " symbol in C/C++ language  !!!
 31             {                            //////////// pow(x,y) -> x^y
 32                 int q = m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
 33                 if(q<m[i][j])
 34                 {
 35                     m[i][j] = q;
 36                     s[i][j] = k;
 37                 }
 38             }
 39         }
 40     }
 41     
 42 }
 43 
 44 //print the multiplication order
 45 void Print_OptimalParens(int s[N+1][N+1],int i,int j)
 46 {
 47     if(i == j)
 48     {
 49         cout<<"A"<<i;    
 50     }
 51     else
 52     {
 53         cout<<"(";
 54         Print_OptimalParens(s,i,s[i][j]);
 55         Print_OptimalParens(s,s[i][j]+1,j);
 56         cout<<")";
 57     }
 58 }
 59 
 60 // To generate matrices, all matrices have "the same" size
 61 void Init_Matrix(long A[][MAX][MAX],int *p)
 62 {    
 63     //read metrix data from file "MetrixData" generated by rand();
 64     ifstream fin("MetrixData.txt");
 65     //vector<int> vec;
 66     long idata;
 67     fin>>idata;
 68     for(int i=1;i<=N;i++)
 69     {
 70         for(int j=0;j<p[i-1];j++) //matrices have different size!
 71         {
 72             for(int k=0;k<p[i];k++)
 73             {                
 74                 A[i][j][k] = idata;
 75                 fin>>idata;    
 76             }
 77         }
 78     }
 79 }
 80 
 81 Mat Matrix_Multiplication(Mat a,Mat b,Mat c)
 82 {
 83     for(int i=0;i<MAX;i++)
 84     {    
 85         for(int j=0;j<MAX;j++)
 86         {
 87             for(int k=0;k<MAX;k++)
 88             {
 89                 c[i][j]+=a[i][k]*b[k][j];
 90             }
 91         }
 92     }
 93     return c;
 94 }
 95 
 96 Mat Matrix_Chain_Multiplication(int i,int j,int s[N+1][N+1])
 97 {
 98     if(i == j)
 99     {
100         return A[i];
101     }
102     else if(i+1 == j)
103     {
104         pp[num] = (Mat)malloc(MAX*MAX*sizeof(long));
105         if(pp[num] == NULL)
106         {
107             cout<<"Can\'t allocate memory!";
108             exit(1);
109         }
110         memset(pp[num],0,MAX*MAX*sizeof(long));
111         return Matrix_Multiplication(A[i],A[j],pp[num ++]);    
112         // num ++; // This line can\'t be excuted !!!
113     }
114     else
115     {    
116         Mat l_Mat = Matrix_Chain_Multiplication(i,s[i][j],s);
117         Mat r_Mat = Matrix_Chain_Multiplication(s[i][j]+1,j,s);
118         
119         pp[num] = (Mat)malloc(MAX*MAX*sizeof(long));
120         if(pp[num] == NULL)
121         {
122             cout<<"Can\'t allocate memory!";
123             exit(1);
124         }
125         memset(pp[num],0,MAX*MAX*sizeof(long));
126         
127         return Matrix_Multiplication(l_Mat,r_Mat,pp[num ++]);
128     }
129 }
130 
131 //print 2_Dim array
132 void Print_Array(long b[][MAX])
133 {
134     for(int i=0;i<MAX;i++)
135     {    
136         for(int j=0;j<MAX;j++)
137         {
138             cout<<b[i][j]<<" ";
139         }
140         cout<<endl;
141     }
142 }
143 
144 void main()
145 {
146     int m[N+1][N+1];
147     int s[N+1][N+1];
148     int p[N+1]={10,8,5,7,6,12,9}; //store the size of matrices
149     Matrix_Chain_Order(p,m,s);
150         
151     // print m[][]
152     for(int i=1;i<=N;i++)
153     {
154         for(int j=i+1;j<=N;j++)
155         {
156             cout<<" m["<<i<<"]["<<j<<"]: "<<m[i][j];
157         }
158         cout<<endl;
159     }
160     //print s[][]
161     for(i=1;i<=N;i++)
162     {
163         for(int j=i+1;j<=N;j++)
164         {
165             cout<<" s["<<i<<"]["<<j<<"]: "<<s[i][j];
166         }
167         cout<<endl;
168     }
169 
170     Print_OptimalParens(s,1,6);
171     cout<<endl;
172 
173     Init_Matrix(A,p);
174     //Print_Array(A[1]);
175     Matrix_Chain_Multiplication(1,6,s);
176     Print_Array(pp[--num]);
177     
178     for(i=0;i<=num;i++)
179     {
180         free(pp[i]);
181     }
182     
183 }

OptimalMatrixMultiplication.cpp