数学建模社区-数学中国

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

---

作者: forcal    时间: 2011-10-24 18:54
标题: 在Lu中创建自定义数据类型，小矩阵乘效率测试
    本例中，我们将自定义矩阵（matrix）类型，基本类型和扩展类型均为matrix（标识矩阵）。
. r( w9 B. U3 B' h3 y5 V8 Q" w% \
    基本要点：

8 d7 _5 e6 @8 x" E" _, @    （1）编写生成矩阵（matrix）的函数NewMatrix和销毁矩阵的函数DelMatrix。

    （2）为自定义类型matrix编写运算符重载函数OpMatrix。

4 g* ]7 {  o* f" I# M# V( K2 Z7 x    （3）用函数LockKey将重载函数OpMatrix注册到Lu，锁定的键的类型即为matrix，要注册为常量，以便于使用。
8 u: D% Q/ S5 g+ J" I: U
    （4）为自定义类型matrix编写其他操作函数（本例未提供）。
; Q. ~! {; v* ^* j
4 a! G' r5 _" u. N    （5）用函数LockKey解锁键matrix（本例中，程序退出时会自动解锁，故可以不用）。

1. #include <windows.h>
   
2. #include <iostream>
   
3. #include <math.h>
   
4. #include "lu32.h"
   # g6 H* Z! M. N3 }. O* O
5. #pragma comment( lib, "lu32.lib" )
   
6. using namespace std;
   
7. //自定义矩阵
     K" Y/ X: K5 k/ B  ~8 Q( a4 b" x' J
8. class myMatrix
   
9. {
   
10. public:
    1 c/ ?+ z0 a+ z5 y$ R9 a: L' P
11.         double *Array;        //数据缓冲区
    
12.         luVOID ArrayLen;        //数据缓冲区长度
    
13.         luVOID Dim[2];        //矩阵维数
    
14.         myMatrix(){Array=NULL; ArrayLen=0; Dim[0]=0; Dim[1]=0;}
    + O7 h  z' ~" g
15.         ~myMatrix()
    - [0 j6 P, V& ^9 k
16.         {
    
17.                 if(Array) delete[] Array;
    " l' D% T- P: @. n  a9 A7 T+ H3 p
18.         }
    # j# z3 \! I! p9 w) U9 ^
19. };
    
20. luKEY Matrix=-1000;        //标识矩阵类型，最终的Matrix由LockKey决定
    6 D$ @0 }& W6 w. |$ z+ q
21. void _stdcall LuMessage(wchar_t *pch)//输出动态库信息，该函数注册到Lu，由Lu二级函数调用
    
22. {
    # H3 {# X4 N, e
23.         wcout<<pch;
    ) _0 x& G; ~$ y+ r* w* `
24. }
    + g0 G& E( O! z* c9 ]( W
25. void _stdcall DelMatrix(void *me)        //用于LockKey函数及InsertKey函数，使Lu能自动销毁myMatrix对象
    
26. {
    
27.         delete (myMatrix *)me;
    3 f- t) W5 S6 J) g3 C3 c8 b
28. }
    
29. myMatrix * _stdcall NewMatrix(luVOID m,luVOID n)        //生成一个myMatrix对象
    $ [  u) \9 k8 y  p% N
30. {
    
31.         myMatrix *pMatrix;
    , b! P/ Y# Z+ P9 b% S# v
32.         luVOID k;
    
33.         double *pa;
    : P: z% T9 ?' S: Z9 o( d2 M# C+ U! W
34.         char keyname[sizeof(luVOID)];
    
35.         void *NowKey;
    
36.         k=m*n;
    1 l, O: J' ^$ W& Q7 K3 O
37.         pMatrix=(myMatrix *)GetBufObj(Matrix,keyname);//先尝试从缓冲区中获取一个矩阵对象
    
38.         if(pMatrix)
    3 Y) `+ @% h7 d9 C  z2 n
39.         {
    
40.                 if(pMatrix->ArrayLen!=k)        //重置矩阵的大小
    7 d5 U+ h; y4 a& L# j2 R4 V* o
41.                 {
    
42.                         pa=new double[k];
    9 i( K/ b& C) K5 l. h
43.                         if(!pa)
    
44.                         {
    
45.                                 DeleteKey(keyname,sizeof(luVOID),Matrix,DelMatrix,1);        //将矩阵对象放回缓冲区
    
46.                                 return NULL;
    ' ^6 s* j' u2 j  v* }
47.                         }
    * W) Z  n/ m2 E4 G* _3 B
48.                         delete[] pMatrix->Array;
    
49.                         pMatrix->Array=pa;
    
50.                 }
    ! o" f1 d& F, H5 l7 [+ h3 Y
51.         }
    / S8 w$ k; @' n
52.         else
    
53.         {
      V8 \% k5 c; |0 W$ X' U- }# Z
54.                 pMatrix=new myMatrix;        //创建矩阵对象
    
55.                 if(!pMatrix) return NULL;
    ( Z7 b8 h! r8 Q# v1 c. P# d* N  a! a
56.                 pMatrix->Array=new double[k];
    
57.                 if(!pMatrix->Array)
    ( ^7 Y+ Q- e" L! Y
58.                 {
    
59.                         delete pMatrix;
    
60.                         return NULL;
    
61.                 }
    
62.                 if(InsertKey((char *)&pMatrix,-1,Matrix,pMatrix,DelMatrix,NULL,0,NowKey))        //将矩阵对象注册到Lu
    
63.                 {
    5 U+ D7 u- C) a3 j  c: ?4 J. z
64.                         delete pMatrix;
    $ G5 p* j3 Q; _! |$ I* K7 T
65.                         return NULL;
    1 W" F& w4 \7 g7 d& a( J
66.                 }
    8 ~- U. _! o) a$ I
67.         }
    % Q1 ^& P- a& P+ P: r4 d
68.         pMatrix->ArrayLen=k; pMatrix->Dim[0]=m; pMatrix->Dim[1]=n;
    
69.         return pMatrix;
    
70. }
    $ R" k6 A- ~* _! O9 M6 a' [
71. LuData _stdcall OpMatrix(luINT mm,LuData *xx,void *hFor,int theOperator)        //运算符重载函数，用于LockKey函数
    $ L5 i+ l5 y, f1 ]
72. {
    
73.         LuData a;
    
74.         myMatrix *pMatrix1,*pMatrix2,*pMatrix3;
    
75.         luVOID i,j,k,m,n,u,v;
    ( Z! }, O: u$ m. x2 H! @
76.         double *pa,*pb,*pc;
    
77.         luMessage pMessage;
      p) ^0 m' r! |6 [9 }% J5 V6 l
78.         wchar_t wchNum[32];
    
79.         char chNum[32];
    
80.         a.BType=luStaData_nil; a.VType=luStaData_nil; a.x=0;
    1 |1 t1 A( ^/ V* E& A" W
81.         switch(theOperator)
    
82.         {
    
83.         case 2:        //重载运算符*
    
84.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
      _2 t7 S: c7 b# m' K5 a
85.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
    
86.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
    
87.                 if(pMatrix1->Dim[1]!=pMatrix2->Dim[0]) break;        //维数不匹配
    ( d6 z0 [% d: ]5 {& ?8 J( v
88.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix2->Dim[1]);        //生成新矩阵
    - {3 f. P  c) M+ [' S6 S
89.                 if(!pMatrix3) break;
    
90.                 pa=pMatrix1->Array; pb=pMatrix2->Array; pc=pMatrix3->Array;
    " R  Q6 T8 J3 |4 u0 j, L
91.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=pMatrix2->Dim[1];
    
92.                 for(i=0; i<m; i++)        //矩阵乘
    
93.                 {
    " h; B9 q5 c; G; m- f  e0 l0 z- w! r
94.                         for(j=0; j<k; j++)
    % I# O3 K1 }1 e& n. o6 g9 V
95.                         {
    
96.                                 u=i*k+j; pc[u]=0.0;
    
97.                                 for (v=0; v<n; v++)
    
98.                                 {
    
99.                                         pc[u]=pc[u]+pa[i*n+v]*pb[v*k+j];
    # N: D' a% O5 _  l$ K
100.                                 }
     
101.                         }
     
102.                 }
     
103.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     2 J& j. o* k1 m7 u9 n' {3 j. x  k
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);
     ' c& z& o- b1 s% r) l+ [
108.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
     $ c! w& f4 A3 u: _
109.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
     * X+ E( s  O* r% l
110.                 if(pMatrix1->Dim[0]!=pMatrix2->Dim[0] || pMatrix1->Dim[1]!=pMatrix2->Dim[1]) break;        //维数不相同
     
111.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix1->Dim[1]);        //生成新矩阵
     ) U" a+ b: J* _/ C+ r+ z* A
112.                 if(!pMatrix3) break;
     & I( R; h$ A! ^: f. l
113.                 for(i=0;i<pMatrix1->ArrayLen;i++) pMatrix3->Array[i]=pMatrix1->Array[i]*pMatrix2->Array[i];        //矩阵点乘
     
114.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
115.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     0 L% P& l2 H  t: P' @9 I0 _( w8 y5 C
116.                 break;
     
117.         case 46:        //重载函数new
     7 e4 O; t8 d3 S7 ]
118.                 if(mm<2) break;
     ! [5 H/ r1 [9 F: N$ F0 ~
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);//生成新矩阵
     
121.                 if(!pMatrix3) break;
     2 D1 G3 H( |% o& d
122.                 for(j=0,i=3;i<=mm;i++,j++)        //赋初值
     + g5 q2 Z* v9 _: D, @6 \' x* u  g0 y
123.                 {
     
124.                         if(j>=pMatrix3->ArrayLen) break;
     
125.                         if((xx+i)->BType!=luStaData_double) break;        //只接受实数参数
     / c0 f5 `* C+ t2 ?7 t6 N# g8 \
126.                         pMatrix3->Array[j]=*(double *)&((xx+i)->x);
     
127.                 }
     
128.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
129.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     $ {. t6 |2 `1 C
130.                 break;
     
131.         case 49:        //重载函数o
     
132.                 pMessage=(luMessage)SearchKey("\0\0\0\0",sizeof(luVOID),luPubKey_User);
     
133.                 if(!pMessage) break;
     4 u* a: f$ r) r7 H7 Q, I* @+ p
134.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     
135.                 if(!pMatrix1) break;        //对象句柄无效，不是矩阵
     ' H0 W# s: [1 Z" O; O8 _, R5 b% x
136.                 pa=pMatrix1->Array;
     4 G& s( ?1 o0 m- M* U
137.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=0;
     
138.                 for(i=0; i<m; i++)        //输出矩阵
     
139.                 {
     2 Z" K: b8 c/ f  j' p1 A
140.                         pMessage(L"\r\n"); k+=2;
     # _% X& k" m0 Z0 z- E! \
141.                         for(j=0; j<n; j++)
     8 g' B0 Q, L) s& Y( r$ d* M5 y
142.                         {
     2 j9 @1 |5 _4 ~  b! [& Z/ l0 p
143.                                 _gcvt_s(chNum,pa[i*n+j],16);
     
144.                                 for(u=0;chNum[u];u++) {wchNum[u]=chNum[u]; k++;}
       w6 i3 k. I3 ?: o! j# a2 M
145.                                 wchNum[u]='\0';
     5 G! F  z6 p3 K5 C& `/ ~
146.                                 pMessage(wchNum); pMessage(L"  "); k+=2;
     ! j  }1 w" Y4 O2 k( L( Z1 K% H
147.                         }
     
148.                 }
     
149.                 pMessage(L"\r\n"); k+=2;
     
150.                 a.BType=luStaData_int64; a.VType=luStaData_int64; a.x=k;        //按函数o的要求，返回输出的字符总数
     9 O- K' u! ~! `9 J
151.                 break;
     
152.         default:
     
153.                 break;
     
154.         }
     
155.         return a;
     
156. }
     3 M! s8 N# |7 m
157. void main(void)
     / \2 F# b9 G6 K: W8 J' Q
158. {
     
159.         void *hFor;                //表达式句柄
     
160.         luINT nPara;                //存放表达式的自变量个数
     
161.         LuData *pPara;                //存放输入自变量的数组指针
     
162.         luINT ErrBegin,ErrEnd;        //表达式编译出错的初始位置和结束位置
     ' N+ J4 R, }2 Z  W9 g
163.         int ErrCode;                //错误代码
     ( g1 U% d6 X& r# A
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.]}";//字符串表达式，矩阵乘
     - T: N' E1 l7 ~" i9 l. _
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;
     3 t& z! h: ?9 I% H6 a" b/ ]
168.         if(!InitLu()) return;        //初始化Lu
     
169.         while(LockKey(Matrix,DelMatrix,OpMatrix)){Matrix--;}        //锁定一个键，用于存储矩阵扩展类型
     
170. 
     ( ~8 i7 {% N+ z) C. L0 o
171.         Val.BType=luStaData_int64; Val.VType=luStaData_int64; Val.x=Matrix;        //定义整数常量
     * {2 Y  T# b& p5 {8 l
172.         SetConst(L"matrix",&Val);        //设置整数常量
     6 U/ d, t1 `9 _4 s0 G* [
173.         InsertKey("\0\0\0\0",4,luPubKey_User,LuMessage,NULL,NULL,1,v); //使Lu运行时可输出函数信息
     0 M! \- [7 j! D1 t
174.         wcout.imbue(locale("chs"));        //设置输出的locale为中文
     
175.   
     
176.         ErrCode=LuCom(ForStr,0,0,0,hFor,nPara,pPara,ErrBegin,ErrEnd); //编译表达式
     
177.         if(ErrCode)
     
178.         {
     : T4 X- n0 A' g
179.                 wcout<<L"表达式有错误！错误代码："<<ErrCode<<endl;
     ) F" ?  O$ r5 v" G0 p
180.         }
     % j, p* G9 y5 \7 n" V9 D* j$ H9 ~
181.         else
       n$ ]+ W) L* i! j7 K$ n2 e  x. i
182.         {
     
183.                 LuCal(hFor,pPara);        //计算表达式的值
     
184.         }
     & w* y) y7 \( O! d7 w- h5 ~, c. u
185.         LockKey(Matrix,NULL,OpMatrix);//解锁键Matrix，本例中，该函数可以不用
     
186.         FreeLu();                        //释放Lu
     
187. }

复制代码

习题：

    （1）自定义矩阵的加、减、左除、右除、点左除等运算，自编测试字符串代码，重新编译运行程序，观察计算结果。
. l: `2 S; p9 k1 b, v
& F8 ^( |$ R4 ]$ ^+ H5 l    （2）小矩阵乘效率测试。编译运行以下Lu字符串代码：

1. main(:a,b,c,d,t,i)=
   : B$ a( o1 u. W/ `2 `# Q# o
2.     a=new[matrix,2,2: 1.,2.,2.,1.],
   , Z" g) A% W2 m) G6 E5 Y# M
3.     b=new[matrix,2,2: 2.,1.,1.,2.],
   3 N+ ^( A% A8 u; j
4.     c=new[matrix,2,2: 2/3.,-1/3.,-1/3.,2/3.],
   
5.     t=clock(),
   
6.     d=a*b, i=0, while{i<1000000, d=d*c*b, i++},
   9 e  `5 }# h- P
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.
   
3. time=0.797 seconds.
   
4. 请按任意键继续. . .

复制代码

Matlab 2009a 代码：

1. a=[1.,2.;2.,1.];
   
2. b=[2.,1.;1.,2.];
   1 A: b) U; @  A5 F8 j/ o
3. c=[2/3.,-1/3.;-1/3.,2/3.];
   
4. tic,
   , M* {7 T. A2 m2 Q+ E9 H! }2 u! q% e
5. d=a*b;
   3 B. q7 N) n9 w2 i2 X4 _0 R
6. for i=1:1000000
   
7.     d=d*c*b;
   
8. end
   # V2 _  i6 }0 e+ D' D( k
9. d,
   
10. toc

复制代码

结果：

1. d =
   
2.      4     5
   
3.      5     4
   
4. Elapsed time is 2.903034 seconds.

复制代码

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

. T  b2 c( W! ?    由以上可以看出，自定义数据类型和系统内置类型有近乎相同的效率。

---

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

Powered by Discuz! X2.5