数学建模社区-数学中国

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

---

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

& B' Q4 W. E; j    基本要点：

    （1）编写生成矩阵（matrix）的函数NewMatrix和销毁矩阵的函数DelMatrix。
+ h$ A5 [7 U+ O3 Y. r5 @8 a
    （2）为自定义类型matrix编写运算符重载函数OpMatrix。
$ ^! B! b, A- C6 b( w
# h  F, {' o" h1 @3 s8 W    （3）用函数LockKey将重载函数OpMatrix注册到Lu，锁定的键的类型即为matrix，要注册为常量，以便于使用。

  k% C+ E* R8 N* I, P0 q7 N1 d    （4）为自定义类型matrix编写其他操作函数（本例未提供）。

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

1. #include <windows.h>
   8 B3 b5 ]% j$ e; J7 W; a) s
2. #include <iostream>
   5 f5 V0 ~: K4 d0 i8 X  ~
3. #include <math.h>
   0 O. \7 w7 a  u* w' ?. o, G% P. |
4. #include "lu32.h"
   ; `5 N0 O1 d( q5 ?
5. #pragma comment( lib, "lu32.lib" )
   
6. using namespace std;
   
7. //自定义矩阵
   7 S3 ?4 ?- Y# K" L# g
8. class myMatrix
   & [3 U* D# l7 S" [- [/ e
9. {
   
10. public:
    " q, C. ^( M; }; ]8 \6 B
11.         double *Array;        //数据缓冲区
    6 \% j/ B( R' I, n) e8 U$ U1 ~# t
12.         luVOID ArrayLen;        //数据缓冲区长度
    ; v2 L  [9 j  G" I
13.         luVOID Dim[2];        //矩阵维数
    
14.         myMatrix(){Array=NULL; ArrayLen=0; Dim[0]=0; Dim[1]=0;}
    ) V7 ]# S( j* j$ H, z' P2 Z
15.         ~myMatrix()
    
16.         {
    7 G0 }+ t8 A+ e9 f- P
17.                 if(Array) delete[] Array;
    
18.         }
    % ]0 t( X5 D' f6 [* S9 h  x
19. };
    
20. luKEY Matrix=-1000;        //标识矩阵类型，最终的Matrix由LockKey决定
    
21. void _stdcall LuMessage(wchar_t *pch)//输出动态库信息，该函数注册到Lu，由Lu二级函数调用
    ( ]) r" j, H5 V* F# R' ]
22. {
    
23.         wcout<<pch;
    8 V- ^1 I) e$ k: H$ X) }7 f# Z% h" l
24. }
    
25. void _stdcall DelMatrix(void *me)        //用于LockKey函数及InsertKey函数，使Lu能自动销毁myMatrix对象
    
26. {
    
27.         delete (myMatrix *)me;
    * g! M% ^$ R- M0 J
28. }
    % E1 ]) y' T) }! C1 f, I. j
29. myMatrix * _stdcall NewMatrix(luVOID m,luVOID n)        //生成一个myMatrix对象
    
30. {
    / F  p* \2 j! n7 L4 s5 v' l( R6 k
31.         myMatrix *pMatrix;
    / @& `  t, }  D
32.         luVOID k;
    3 F% S2 g3 F" P- @
33.         double *pa;
    : @- T# r0 E8 x$ m( Y6 l+ w3 g
34.         char keyname[sizeof(luVOID)];
    & p( Y7 P/ m0 Z# h0 S. g2 i  w
35.         void *NowKey;
    
36.         k=m*n;
    
37.         pMatrix=(myMatrix *)GetBufObj(Matrix,keyname);//先尝试从缓冲区中获取一个矩阵对象
    
38.         if(pMatrix)
    ) k$ P0 w9 a- s1 n# ~; o. D
39.         {
    
40.                 if(pMatrix->ArrayLen!=k)        //重置矩阵的大小
    
41.                 {
    
42.                         pa=new double[k];
    
43.                         if(!pa)
    
44.                         {
    
45.                                 DeleteKey(keyname,sizeof(luVOID),Matrix,DelMatrix,1);        //将矩阵对象放回缓冲区
    
46.                                 return NULL;
    
47.                         }
    
48.                         delete[] pMatrix->Array;
    
49.                         pMatrix->Array=pa;
    
50.                 }
    9 T8 P7 ~& |- w1 T# x
51.         }
    
52.         else
    ) G; |0 t3 H$ p% ^
53.         {
    
54.                 pMatrix=new myMatrix;        //创建矩阵对象
    " E) i2 O( N8 D8 }" @
55.                 if(!pMatrix) return NULL;
    
56.                 pMatrix->Array=new double[k];
    0 v: _5 x7 d- H: K
57.                 if(!pMatrix->Array)
    
58.                 {
    $ s8 E1 B7 B" z4 O, C  O% a
59.                         delete pMatrix;
    
60.                         return NULL;
    7 F* G% Y1 x; }1 U
61.                 }
    
62.                 if(InsertKey((char *)&pMatrix,-1,Matrix,pMatrix,DelMatrix,NULL,0,NowKey))        //将矩阵对象注册到Lu
    $ |9 b( |# T: X1 `& j- G) M6 [2 j
63.                 {
    & q0 p4 I9 ]. @; z8 i
64.                         delete pMatrix;
    9 U1 v0 J% l$ a5 ~" y/ T
65.                         return NULL;
    
66.                 }
    1 a* a% |5 {* N$ g& Q& q) G
67.         }
    3 W$ Y2 Y& ]0 d, M; W: ?
68.         pMatrix->ArrayLen=k; pMatrix->Dim[0]=m; pMatrix->Dim[1]=n;
    % j2 N' o, x6 n$ T5 Z- T2 ^
69.         return pMatrix;
    * o4 j: c( c. H% v  Y6 U6 x* T
70. }
    ' x# Y" S% {* |
71. LuData _stdcall OpMatrix(luINT mm,LuData *xx,void *hFor,int theOperator)        //运算符重载函数，用于LockKey函数
    
72. {
    
73.         LuData a;
    
74.         myMatrix *pMatrix1,*pMatrix2,*pMatrix3;
    ; m$ g3 y9 |8 J
75.         luVOID i,j,k,m,n,u,v;
      r/ F# @( T3 n( K0 H* q8 ^+ F
76.         double *pa,*pb,*pc;
    ) I/ o9 x! _% f2 k/ y' s# a
77.         luMessage pMessage;
    # Z: d, c& S: I: r" h
78.         wchar_t wchNum[32];
    
79.         char chNum[32];
    : E" P# |2 Q- x  O
80.         a.BType=luStaData_nil; a.VType=luStaData_nil; a.x=0;
      `- x7 {: c7 z# }) ~& s
81.         switch(theOperator)
    / H# B$ _1 B% g2 v6 M
82.         {
    / W1 E* b* h, U
83.         case 2:        //重载运算符*
    
84.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
    : b7 g# f$ Q/ ~5 Z  Z
85.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
    
86.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
    9 p- _$ w# Y) \; R1 i( p5 a6 d
87.                 if(pMatrix1->Dim[1]!=pMatrix2->Dim[0]) break;        //维数不匹配
    
88.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix2->Dim[1]);        //生成新矩阵
    
89.                 if(!pMatrix3) break;
    9 Z3 R* B! R( ^9 _& U, S: G
90.                 pa=pMatrix1->Array; pb=pMatrix2->Array; pc=pMatrix3->Array;
    9 o) w# G' y, N' r9 ?
91.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=pMatrix2->Dim[1];
    4 p- v8 s, K) u2 k* c& x# ]8 c4 N5 S
92.                 for(i=0; i<m; i++)        //矩阵乘
    & _. \  f2 P3 K- @- h' t0 K8 L$ t' a
93.                 {
    0 l$ x9 |4 B; b; O$ d" `# C9 R0 j5 t0 \2 P  I
94.                         for(j=0; j<k; j++)
    
95.                         {
    9 @) ?+ C* ^7 z) Z) s
96.                                 u=i*k+j; pc[u]=0.0;
    
97.                                 for (v=0; v<n; v++)
    ' R# i% [" a: u# r5 H
98.                                 {
    
99.                                         pc[u]=pc[u]+pa[i*n+v]*pb[v*k+j];
    . s% t0 t+ l. Q$ @& U# p
100.                                 }
     
101.                         }
     + l9 H2 `5 n# I: V1 t' U; `( I
102.                 }
     # N3 ~- L4 n# f- _
103.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     6 C, |. o- g3 x% P+ z4 c
104.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     3 B" c: c; ^) Z, {& j* }
105.                 break;
     0 G6 b$ k* m, I4 m; @! h; k
106.         case 25:        //重载运算符.*
     ! U. W8 W" k0 X0 W
107.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     3 Y7 q; L% n* w; p2 L) l& X8 i+ l
108.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
     
109.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
     
110.                 if(pMatrix1->Dim[0]!=pMatrix2->Dim[0] || pMatrix1->Dim[1]!=pMatrix2->Dim[1]) break;        //维数不相同
     " [3 K* f9 s3 m' A) E
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];        //矩阵点乘
     
114.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
115.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     / p/ S9 G# _( a9 X- n
116.                 break;
     
117.         case 46:        //重载函数new
     
118.                 if(mm<2) break;
     " `1 O+ `3 _2 ^1 U5 N$ c# @% K2 z  E8 K/ @
119.                 if((xx+1)->x<1 || (xx+2)->x<1 || (xx+1)->BType!=luStaData_int64 || (xx+2)->BType!=luStaData_int64) break;
     2 C& T! h+ R" ?
120.                 pMatrix3=NewMatrix((luVOID)(xx+1)->x,(luVOID)(xx+2)->x);//生成新矩阵
     7 Q) F" `# i  a5 z' V' ?! s
121.                 if(!pMatrix3) break;
     $ X$ H$ u; i# d( Y. E. x0 P' H/ f
122.                 for(j=0,i=3;i<=mm;i++,j++)        //赋初值
     % H( e" ~# F% V  w; f; x- k
123.                 {
     
124.                         if(j>=pMatrix3->ArrayLen) break;
     
125.                         if((xx+i)->BType!=luStaData_double) break;        //只接受实数参数
     
126.                         pMatrix3->Array[j]=*(double *)&((xx+i)->x);
     7 x  w' c3 X- o1 @* g, k9 K
127.                 }
     , d' t+ i( ?7 e5 T- N* `" d& ?
128.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
129.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     ' J' {0 x% r/ B
130.                 break;
     % a8 K) |' H6 ~6 v4 U, g& z0 @! {
131.         case 49:        //重载函数o
     
132.                 pMessage=(luMessage)SearchKey("\0\0\0\0",sizeof(luVOID),luPubKey_User);
     
133.                 if(!pMessage) break;
     
134.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     , m5 s, C3 T3 P5 T0 F3 D2 T
135.                 if(!pMatrix1) break;        //对象句柄无效，不是矩阵
     
136.                 pa=pMatrix1->Array;
     # ~+ E! X. m* ?
137.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=0;
     
138.                 for(i=0; i<m; i++)        //输出矩阵
     : T. G& T) N5 T1 M
139.                 {
     6 p3 ]. b0 _1 Q& T/ b7 J8 M
140.                         pMessage(L"\r\n"); k+=2;
     3 Q( {0 ]0 Y0 }5 o9 H$ P
141.                         for(j=0; j<n; j++)
     
142.                         {
     
143.                                 _gcvt_s(chNum,pa[i*n+j],16);
     
144.                                 for(u=0;chNum[u];u++) {wchNum[u]=chNum[u]; k++;}
     2 E& i: J* [% Z: P* Q
145.                                 wchNum[u]='\0';
     
146.                                 pMessage(wchNum); pMessage(L"  "); k+=2;
     ' a! I  A# U$ |" ~
147.                         }
     
148.                 }
     
149.                 pMessage(L"\r\n"); k+=2;
     
150.                 a.BType=luStaData_int64; a.VType=luStaData_int64; a.x=k;        //按函数o的要求，返回输出的字符总数
     - C" j2 m, n: ]
151.                 break;
     
152.         default:
     # J4 e. q( P% I- N6 @
153.                 break;
     
154.         }
     0 E; }2 t! l+ `5 o; |+ `# D- l
155.         return a;
     6 w* E- A; ^/ o3 _$ U9 G' y
156. }
     
157. void main(void)
     ! e  H  n7 \# q& f, W/ r
158. {
     & }/ j! z  x# z: X% h7 k. t/ ?
159.         void *hFor;                //表达式句柄
     
160.         luINT nPara;                //存放表达式的自变量个数
     
161.         LuData *pPara;                //存放输入自变量的数组指针
     
162.         luINT ErrBegin,ErrEnd;        //表达式编译出错的初始位置和结束位置
     
163.         int ErrCode;                //错误代码
     9 k) m1 ~5 z0 }7 M
164.         void *v;
     $ q9 u" V( |3 d! L
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.]}";//字符串表达式，矩阵乘
     / D  O) f$ ?4 d2 i. ^/ n1 d2 H# t1 k1 K
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;
     
168.         if(!InitLu()) return;        //初始化Lu
     % e4 v& g, t1 x5 u. S' z
169.         while(LockKey(Matrix,DelMatrix,OpMatrix)){Matrix--;}        //锁定一个键，用于存储矩阵扩展类型
     3 x  J3 W9 B* v6 E4 W
170. 
     7 J2 t% H. ?* `3 b% T) B
171.         Val.BType=luStaData_int64; Val.VType=luStaData_int64; Val.x=Matrix;        //定义整数常量
     # c& J& k; t3 U) L6 b* h9 t6 a
172.         SetConst(L"matrix",&Val);        //设置整数常量
     7 y/ I$ q7 l' P; k& r
173.         InsertKey("\0\0\0\0",4,luPubKey_User,LuMessage,NULL,NULL,1,v); //使Lu运行时可输出函数信息
     
174.         wcout.imbue(locale("chs"));        //设置输出的locale为中文
     
175.   
     
176.         ErrCode=LuCom(ForStr,0,0,0,hFor,nPara,pPara,ErrBegin,ErrEnd); //编译表达式
     
177.         if(ErrCode)
     
178.         {
     " J; i8 ?/ J; M- \1 f
179.                 wcout<<L"表达式有错误！错误代码："<<ErrCode<<endl;
     ) M4 v3 \2 T" e% n1 m
180.         }
     
181.         else
     
182.         {
     / _0 p) l- _8 S, g' g
183.                 LuCal(hFor,pPara);        //计算表达式的值
     , b7 l- u2 ~7 L6 j  Q
184.         }
     
185.         LockKey(Matrix,NULL,OpMatrix);//解锁键Matrix，本例中，该函数可以不用
     
186.         FreeLu();                        //释放Lu
     
187. }

复制代码

习题：
- P4 i; J7 `9 p, Y
    （1）自定义矩阵的加、减、左除、右除、点左除等运算，自编测试字符串代码，重新编译运行程序，观察计算结果。
% v# l5 t$ ^9 _; \$ m  I
$ `% B9 F; _& N$ g    （2）小矩阵乘效率测试。编译运行以下Lu字符串代码：

1. main(:a,b,c,d,t,i)=
   5 b$ `8 |. B) J) o1 @3 p; R
2.     a=new[matrix,2,2: 1.,2.,2.,1.],
   " _7 w6 U+ J' P& N- G3 x2 D
3.     b=new[matrix,2,2: 2.,1.,1.,2.],
   6 l- r5 s! b. c% L0 F
4.     c=new[matrix,2,2: 2/3.,-1/3.,-1/3.,2/3.],
   : o  i4 k& B- S
5.     t=clock(),
   
6.     d=a*b, i=0, while{i<1000000, d=d*c*b, i++},
   & U% \& o9 B3 r3 ~6 Z4 s
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.
   6 b$ n8 v4 `* z0 n
2. 5. 4.
   
3. time=0.797 seconds.
   . ?' I* q) G8 v( q
4. 请按任意键继续. . .

复制代码

Matlab 2009a 代码：

1. a=[1.,2.;2.,1.];
   + l4 n9 i: r9 J5 y4 ]
2. b=[2.,1.;1.,2.];
   
3. c=[2/3.,-1/3.;-1/3.,2/3.];
   - Z6 }5 A; l& c6 S% B+ m
4. tic,
   % _2 P1 S( ]  P/ l
5. d=a*b;
   
6. for i=1:1000000
   * F4 j- L; q5 {5 b
7.     d=d*c*b;
   # |7 Y5 ]% v% a
8. end
   & g5 Z7 g) o: ~$ ^* k/ p
9. d,
   " @1 n! g, y7 ?$ _% q1 n4 _) T
10. toc

复制代码

结果：

1. d =
   
2.      4     5
   ( g; |! r1 c. ?5 T4 ]
3.      5     4
   
4. Elapsed time is 2.903034 seconds.

复制代码

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

7 b; q  y& h8 T8 m- K+ q' H    由以上可以看出，自定义数据类型和系统内置类型有近乎相同的效率。

---

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

Powered by Discuz! X2.5