数学建模社区-数学中国

标题: 在Lu中创建自定义数据类型，小矩阵乘效率测试 [打印本页]

---

作者: forcal    时间: 2011-10-24 18:54
标题: 在Lu中创建自定义数据类型，小矩阵乘效率测试
    本例中，我们将自定义矩阵（matrix）类型，基本类型和扩展类型均为matrix（标识矩阵）。

, Z1 M: s- y3 L( i4 C6 L. S+ U    基本要点：
  o* V( p. y" [0 ]
( y2 K3 u3 `8 i" v! d  y! [    （1）编写生成矩阵（matrix）的函数NewMatrix和销毁矩阵的函数DelMatrix。

8 d$ e/ b( `3 F: @    （2）为自定义类型matrix编写运算符重载函数OpMatrix。

) i. p9 a* {2 F8 E0 q    （3）用函数LockKey将重载函数OpMatrix注册到Lu，锁定的键的类型即为matrix，要注册为常量，以便于使用。
) M5 ]" Y9 ]" F8 K/ W  k1 L% q8 ^
    （4）为自定义类型matrix编写其他操作函数（本例未提供）。

    （5）用函数LockKey解锁键matrix（本例中，程序退出时会自动解锁，故可以不用）。

1. #include <windows.h>
   ( a* _# ?) q2 q
2. #include <iostream>
   + M+ p  w# \9 v$ P: ^1 J6 v1 }9 D
3. #include <math.h>
   , P7 i! y& f' U6 ^; w. S
4. #include "lu32.h"
   3 F  `) X( v$ u5 @' S' X9 o" |
5. #pragma comment( lib, "lu32.lib" )
   
6. using namespace std;
   
7. //自定义矩阵
   - `- N# @; x; O: w
8. class myMatrix
   
9. {
   0 w" W9 A. e1 b* G. N+ Y! H3 U3 F7 y! G1 X
10. public:
    7 Q" S+ P; m5 k' W. j/ w
11.         double *Array;        //数据缓冲区
      k* k4 h; B1 `; k) l4 l, g
12.         luVOID ArrayLen;        //数据缓冲区长度
    
13.         luVOID Dim[2];        //矩阵维数
    + v0 T2 a/ f6 x+ K$ ]& ~
14.         myMatrix(){Array=NULL; ArrayLen=0; Dim[0]=0; Dim[1]=0;}
    
15.         ~myMatrix()
    2 g. x( k8 T4 r, i* n% `
16.         {
    
17.                 if(Array) delete[] Array;
    
18.         }
    
19. };
    
20. luKEY Matrix=-1000;        //标识矩阵类型，最终的Matrix由LockKey决定
    
21. void _stdcall LuMessage(wchar_t *pch)//输出动态库信息，该函数注册到Lu，由Lu二级函数调用
    
22. {
    ! k) k3 Q' ?! S% u! v; f1 |
23.         wcout<<pch;
    
24. }
    
25. void _stdcall DelMatrix(void *me)        //用于LockKey函数及InsertKey函数，使Lu能自动销毁myMatrix对象
    0 i) K6 k' u) S) _. Z6 m
26. {
    
27.         delete (myMatrix *)me;
    
28. }
    
29. myMatrix * _stdcall NewMatrix(luVOID m,luVOID n)        //生成一个myMatrix对象
    
30. {
    % }0 f2 ^6 ]" M
31.         myMatrix *pMatrix;
    7 ]: V2 o. _9 y
32.         luVOID k;
    
33.         double *pa;
    : }4 E1 ^. G$ q2 T6 V+ L2 S! e* [; h& D
34.         char keyname[sizeof(luVOID)];
    
35.         void *NowKey;
      q, l! e& a1 D' {; I2 t5 ]0 x7 }
36.         k=m*n;
    
37.         pMatrix=(myMatrix *)GetBufObj(Matrix,keyname);//先尝试从缓冲区中获取一个矩阵对象
    2 s/ ?6 {7 K! ^, N
38.         if(pMatrix)
    
39.         {
    : G+ E# @. M4 Q9 {7 N' T, q' V
40.                 if(pMatrix->ArrayLen!=k)        //重置矩阵的大小
    3 S2 p1 I& z5 G* k/ z9 u( b8 t
41.                 {
    
42.                         pa=new double[k];
    $ i" S: h$ z0 k6 T+ d1 t# p
43.                         if(!pa)
    ) C4 R7 y" @4 M$ l% q7 O/ o
44.                         {
    4 Q) R: j1 S' H7 k4 w( e7 @
45.                                 DeleteKey(keyname,sizeof(luVOID),Matrix,DelMatrix,1);        //将矩阵对象放回缓冲区
    
46.                                 return NULL;
    
47.                         }
    . @4 N$ k3 }  q9 p: w
48.                         delete[] pMatrix->Array;
    * g5 v* {0 G' S+ R# Z$ Y5 ]1 V6 t0 s
49.                         pMatrix->Array=pa;
    4 K9 R3 x/ d+ Z  o9 a! V
50.                 }
    ) `! h3 @1 t" M; g
51.         }
    
52.         else
    
53.         {
    ) g/ T0 @) F, O. S) D8 \$ r
54.                 pMatrix=new myMatrix;        //创建矩阵对象
    
55.                 if(!pMatrix) return NULL;
    
56.                 pMatrix->Array=new double[k];
    
57.                 if(!pMatrix->Array)
    
58.                 {
    6 ^$ F( u* V9 O# [" d
59.                         delete pMatrix;
    
60.                         return NULL;
    
61.                 }
    
62.                 if(InsertKey((char *)&pMatrix,-1,Matrix,pMatrix,DelMatrix,NULL,0,NowKey))        //将矩阵对象注册到Lu
    ' x# t, X3 R0 e# Y
63.                 {
    
64.                         delete pMatrix;
    
65.                         return NULL;
    
66.                 }
    " A, J; [/ `. x4 y& T+ R# n' `7 f
67.         }
    
68.         pMatrix->ArrayLen=k; pMatrix->Dim[0]=m; pMatrix->Dim[1]=n;
    ) s' k: O! R" J" {! z7 I7 `, H
69.         return pMatrix;
    * L6 q/ E- \, j# p, m- R
70. }
    
71. LuData _stdcall OpMatrix(luINT mm,LuData *xx,void *hFor,int theOperator)        //运算符重载函数，用于LockKey函数
    0 o5 b0 W2 T" k8 i( U9 I
72. {
    . g- i& ^% t8 g+ S) {! j
73.         LuData a;
    2 q6 f9 E2 ~+ Z
74.         myMatrix *pMatrix1,*pMatrix2,*pMatrix3;
    
75.         luVOID i,j,k,m,n,u,v;
    
76.         double *pa,*pb,*pc;
    ; A2 |6 P& B. P0 [+ z- u5 j
77.         luMessage pMessage;
    6 T  q' Z2 m- @/ ?& _9 x! e- v
78.         wchar_t wchNum[32];
    : j+ V4 C) N+ O% A6 K9 A, A" S/ Q
79.         char chNum[32];
    
80.         a.BType=luStaData_nil; a.VType=luStaData_nil; a.x=0;
    9 V% a. ~' G2 I' Y' s
81.         switch(theOperator)
    ) \# \/ L* ~) s" r
82.         {
    ' m) _( _  b. u/ u; h
83.         case 2:        //重载运算符*
    . H1 J1 b5 R# K: L, d/ ?* _
84.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
    
85.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
    
86.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
    ! p: ]; e/ Y1 D% h, W
87.                 if(pMatrix1->Dim[1]!=pMatrix2->Dim[0]) break;        //维数不匹配
    
88.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix2->Dim[1]);        //生成新矩阵
    ! K& b; l  I6 I
89.                 if(!pMatrix3) break;
    
90.                 pa=pMatrix1->Array; pb=pMatrix2->Array; pc=pMatrix3->Array;
    ! }: D8 |8 o% L0 g9 W) n
91.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=pMatrix2->Dim[1];
    - H" X  x( p2 a# A7 g; N
92.                 for(i=0; i<m; i++)        //矩阵乘
    7 Z, x! A9 \5 t" Y0 E8 C$ K
93.                 {
    
94.                         for(j=0; j<k; j++)
    
95.                         {
    % J4 L* F, N& [4 U6 X
96.                                 u=i*k+j; pc[u]=0.0;
    
97.                                 for (v=0; v<n; v++)
    
98.                                 {
    / S, u" [! ~% {( {. t
99.                                         pc[u]=pc[u]+pa[i*n+v]*pb[v*k+j];
    0 J* S) f- @8 _, }. h. X
100.                                 }
     ) W9 ]  w% X, W: T. A
101.                         }
     
102.                 }
     1 x, g2 @) d6 i
103.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
104.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     
105.                 break;
     
106.         case 25:        //重载运算符.*
     
107.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     
108.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
     
109.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
     , h3 z# n7 v5 L) }3 I& W$ m
110.                 if(pMatrix1->Dim[0]!=pMatrix2->Dim[0] || pMatrix1->Dim[1]!=pMatrix2->Dim[1]) break;        //维数不相同
     5 ~& s6 m0 P# }: \
111.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix1->Dim[1]);        //生成新矩阵
     
112.                 if(!pMatrix3) break;
     
113.                 for(i=0;i<pMatrix1->ArrayLen;i++) pMatrix3->Array[i]=pMatrix1->Array[i]*pMatrix2->Array[i];        //矩阵点乘
     $ y8 l7 B# p4 k. G, s8 H
114.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
115.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     8 t9 z" C3 J7 O2 ?3 T: F
116.                 break;
     * {5 q% V9 k4 R  Q- r
117.         case 46:        //重载函数new
     
118.                 if(mm<2) break;
     
119.                 if((xx+1)->x<1 || (xx+2)->x<1 || (xx+1)->BType!=luStaData_int64 || (xx+2)->BType!=luStaData_int64) break;
     
120.                 pMatrix3=NewMatrix((luVOID)(xx+1)->x,(luVOID)(xx+2)->x);//生成新矩阵
       F5 N( x3 H: `" v/ m: x9 s
121.                 if(!pMatrix3) break;
     
122.                 for(j=0,i=3;i<=mm;i++,j++)        //赋初值
     " x0 i, a) L+ b) x/ T4 x+ _/ D
123.                 {
     
124.                         if(j>=pMatrix3->ArrayLen) break;
     
125.                         if((xx+i)->BType!=luStaData_double) break;        //只接受实数参数
     
126.                         pMatrix3->Array[j]=*(double *)&((xx+i)->x);
     
127.                 }
     
128.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     1 A5 F( s+ Q' H9 r: ~. ^" g
129.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     
130.                 break;
     
131.         case 49:        //重载函数o
     - k+ |, @% |" l9 u. @" x
132.                 pMessage=(luMessage)SearchKey("\0\0\0\0",sizeof(luVOID),luPubKey_User);
     ' r0 l  W* l; N1 Y% o% Q0 C% l2 l' |) ?
133.                 if(!pMessage) break;
     
134.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     
135.                 if(!pMatrix1) break;        //对象句柄无效，不是矩阵
       \: s( D2 g: S+ F( d
136.                 pa=pMatrix1->Array;
     
137.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=0;
     ) j/ {& Z- G' c0 {! c6 c
138.                 for(i=0; i<m; i++)        //输出矩阵
     
139.                 {
     ( A: e6 [3 N0 @/ E6 I
140.                         pMessage(L"\r\n"); k+=2;
     
141.                         for(j=0; j<n; j++)
     / a& |7 Y; b6 q3 B  P$ G9 }
142.                         {
     ' Q3 s& F: r  C8 F
143.                                 _gcvt_s(chNum,pa[i*n+j],16);
     
144.                                 for(u=0;chNum[u];u++) {wchNum[u]=chNum[u]; k++;}
       e; K2 Q9 ^1 b) s& S
145.                                 wchNum[u]='\0';
     6 ?$ n5 d% X! ?8 I1 H2 ~
146.                                 pMessage(wchNum); pMessage(L"  "); k+=2;
     9 y& F. m) c% v" l' _
147.                         }
     - U, P+ C$ p- o+ @' |1 p
148.                 }
     
149.                 pMessage(L"\r\n"); k+=2;
     
150.                 a.BType=luStaData_int64; a.VType=luStaData_int64; a.x=k;        //按函数o的要求，返回输出的字符总数
     
151.                 break;
     8 Q8 H. b* I$ D7 A2 A
152.         default:
     
153.                 break;
     
154.         }
     % n! {" M; o; f/ s; Z
155.         return a;
     
156. }
     7 m. A' j% K; k( v, z0 W
157. void main(void)
     
158. {
     3 e6 l+ h7 v' v% k+ v. s" b
159.         void *hFor;                //表达式句柄
     0 I0 `* w$ z+ Y7 |6 l2 a
160.         luINT nPara;                //存放表达式的自变量个数
     4 d# X% m2 E; O: h9 a( C6 |6 Y( N
161.         LuData *pPara;                //存放输入自变量的数组指针
     
162.         luINT ErrBegin,ErrEnd;        //表达式编译出错的初始位置和结束位置
     
163.         int ErrCode;                //错误代码
     
164.         void *v;
     
165.         wchar_t ForStr[]=L"o{new[matrix,2,3: 0.,1.,2.;3.,4.,5.]*new[matrix,3,2: 1.,2.;3.,4.;5.,6.]}";//字符串表达式，矩阵乘
     
166.         //wchar_t ForStr[]=L"o{new[matrix,2,3: 0.,1.,2.;3.,4.,5.].*new[matrix,2,3: 1.,2.,3.;4.,5.,6.]}";//字符串表达式，矩阵点乘
     
167.         LuData Val;
     $ m0 e# j8 C& K/ M( I' M
168.         if(!InitLu()) return;        //初始化Lu
     
169.         while(LockKey(Matrix,DelMatrix,OpMatrix)){Matrix--;}        //锁定一个键，用于存储矩阵扩展类型
     
170. 
     , Y9 F' e$ r# n) d. Z# i' y
171.         Val.BType=luStaData_int64; Val.VType=luStaData_int64; Val.x=Matrix;        //定义整数常量
     ; x/ t9 w# W  O$ R4 b# H
172.         SetConst(L"matrix",&Val);        //设置整数常量
     - C" v$ t8 D' b- g" B5 r
173.         InsertKey("\0\0\0\0",4,luPubKey_User,LuMessage,NULL,NULL,1,v); //使Lu运行时可输出函数信息
     
174.         wcout.imbue(locale("chs"));        //设置输出的locale为中文
     6 v/ p6 _& Q$ N3 Y, Y& d/ L3 z
175.   
     ' ?% {/ a' f+ I5 G" R
176.         ErrCode=LuCom(ForStr,0,0,0,hFor,nPara,pPara,ErrBegin,ErrEnd); //编译表达式
     / c4 c/ W6 w- r) f2 h0 Y- }
177.         if(ErrCode)
     
178.         {
     
179.                 wcout<<L"表达式有错误！错误代码："<<ErrCode<<endl;
     1 m' D* M( e5 P/ l- ]
180.         }
     
181.         else
     ( B2 j) s% G8 P7 A% j
182.         {
     & H8 S/ `( }$ J, G
183.                 LuCal(hFor,pPara);        //计算表达式的值
     6 ^, C( d) z8 S, Q, P
184.         }
     
185.         LockKey(Matrix,NULL,OpMatrix);//解锁键Matrix，本例中，该函数可以不用
     2 w% A( r% Q: k8 K' E
186.         FreeLu();                        //释放Lu
     5 h$ \2 b" Z6 s  X3 ?
187. }

复制代码

习题：

    （1）自定义矩阵的加、减、左除、右除、点左除等运算，自编测试字符串代码，重新编译运行程序，观察计算结果。
8 f- K  H7 K8 u( A0 U/ Z9 _% Z
0 J: v. l6 i1 W5 `9 f" v! Y: q4 b    （2）小矩阵乘效率测试。编译运行以下Lu字符串代码：

1. main(:a,b,c,d,t,i)=
   9 J& u4 e! ?% ]" C3 t' K
2.     a=new[matrix,2,2: 1.,2.,2.,1.],
   
3.     b=new[matrix,2,2: 2.,1.,1.,2.],
   . W- b0 y% h$ ^1 |
4.     c=new[matrix,2,2: 2/3.,-1/3.,-1/3.,2/3.],
   
5.     t=clock(),
   1 h) j4 G$ U3 H% Z0 k
6.     d=a*b, i=0, while{i<1000000, d=d*c*b, i++},
   . h  O  w, B4 k8 O" Q! E$ C! `9 W, j: U
7.     o{d, "time=",[clock()-t]/1000.," seconds.\r\n"}

复制代码

C/C++中的字符串定义为：

1. wchar_t ForStr[]=L"main(:a,b,c,d,t,i)= a=new[matrix,2,2: 1.,2.,2.,1.], b=new[matrix,2,2: 2.,1.,1.,2.], c=new[matrix,2,2: 2/3.,-1/3.,-1/3.,2/3.], t=clock(), d=a*b, i=0, while{i<1000000, d=d*c*b, i++}, o{d, \"time=\",[clock()-t]/1000.,\" seconds.\r\n\"}";//字符串表达式

复制代码

结果：

1. 4. 5.
   
2. 5. 4.
   1 i* _7 r! I- ?2 [
3. time=0.797 seconds.
   
4. 请按任意键继续. . .

复制代码

Matlab 2009a 代码：

1. a=[1.,2.;2.,1.];
   
2. b=[2.,1.;1.,2.];
   
3. c=[2/3.,-1/3.;-1/3.,2/3.];
   ; C- B# l6 d& n* K/ }* o4 {
4. tic,
   
5. d=a*b;
   . s5 C5 B. ^* @4 ~! H$ X" k
6. for i=1:1000000
   
7.     d=d*c*b;
   . C) ^2 f. w9 C2 u
8. end
   
9. d,
   % s9 b* g" M6 ~  w7 a
10. toc

复制代码

结果：

1. d =
   
2.      4     5
   
3.      5     4
   " O: v1 E6 f7 w# A) ]
4. Elapsed time is 2.903034 seconds.

复制代码

本例矩阵乘效率测试，Lu的速度超过了Matlab，主要在于Lu有更高的动态对象管理效率。

    由以上可以看出，自定义数据类型和系统内置类型有近乎相同的效率。

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5