数学建模社区-数学中国

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

---

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

' Z* ?# a+ G' l: L3 z    基本要点：
9 U) F9 j# Z4 B$ u
% U. M  F% r" N( ]+ L    （1）编写生成矩阵（matrix）的函数NewMatrix和销毁矩阵的函数DelMatrix。

    （2）为自定义类型matrix编写运算符重载函数OpMatrix。
4 l2 Q* i& o# \# g# K  X- E
+ `( }4 E& R* m& {$ ?7 R    （3）用函数LockKey将重载函数OpMatrix注册到Lu，锁定的键的类型即为matrix，要注册为常量，以便于使用。

    （4）为自定义类型matrix编写其他操作函数（本例未提供）。

/ V4 ]# l+ D% M% v& y0 R1 t    （5）用函数LockKey解锁键matrix（本例中，程序退出时会自动解锁，故可以不用）。

1. #include <windows.h>
   + B$ K4 Z: k0 F  |
2. #include <iostream>
   / ~$ p$ ?2 P8 r% E( a3 v2 i3 J5 k
3. #include <math.h>
   
4. #include "lu32.h"
   ( p( o' D; U6 X; s
5. #pragma comment( lib, "lu32.lib" )
   
6. using namespace std;
   
7. //自定义矩阵
   & n! H7 c, y) N( |5 N
8. class myMatrix
   
9. {
   , D( H: C7 W/ m0 H: F# ?
10. public:
    ' @2 b% d, Z+ l% B% o' C
11.         double *Array;        //数据缓冲区
    ! u) W! u9 j' P) s$ R# z1 t. [8 J
12.         luVOID ArrayLen;        //数据缓冲区长度
    0 E7 f5 b- }% \* [, k
13.         luVOID Dim[2];        //矩阵维数
    
14.         myMatrix(){Array=NULL; ArrayLen=0; Dim[0]=0; Dim[1]=0;}
    & }5 E4 Q; M' R" M  z7 F9 C
15.         ~myMatrix()
    : G9 X, Q) f% e+ v- [3 a) S
16.         {
    & N0 j( @$ o, H  i
17.                 if(Array) delete[] Array;
    
18.         }
    % S3 k* y+ u4 E1 S8 y$ I
19. };
    
20. luKEY Matrix=-1000;        //标识矩阵类型，最终的Matrix由LockKey决定
    1 |' R6 h2 x( F& x0 a6 E3 m
21. void _stdcall LuMessage(wchar_t *pch)//输出动态库信息，该函数注册到Lu，由Lu二级函数调用
    
22. {
    & {0 z) b/ B) f- V5 t$ b4 R  ^
23.         wcout<<pch;
    " b/ [% u  u0 A6 X: h; [: w$ p/ g
24. }
    ( v5 t8 N" D; x
25. void _stdcall DelMatrix(void *me)        //用于LockKey函数及InsertKey函数，使Lu能自动销毁myMatrix对象
    ' M  u2 Q  k# R% t+ C& X
26. {
    
27.         delete (myMatrix *)me;
    3 v8 e* O- \. A3 u
28. }
    2 Q, a4 h9 O* a" f- h  n: i$ p
29. myMatrix * _stdcall NewMatrix(luVOID m,luVOID n)        //生成一个myMatrix对象
      z- B) N  y/ E" c
30. {
    1 K/ m2 |/ m' I
31.         myMatrix *pMatrix;
      ?; G# N0 h+ q8 M
32.         luVOID k;
    0 B3 |- B7 W& O& U
33.         double *pa;
    
34.         char keyname[sizeof(luVOID)];
    
35.         void *NowKey;
    
36.         k=m*n;
    
37.         pMatrix=(myMatrix *)GetBufObj(Matrix,keyname);//先尝试从缓冲区中获取一个矩阵对象
    
38.         if(pMatrix)
    ' w6 \% d0 Q+ o) P" g' X
39.         {
    
40.                 if(pMatrix->ArrayLen!=k)        //重置矩阵的大小
    6 G: k4 F6 v9 F! r
41.                 {
    
42.                         pa=new double[k];
    
43.                         if(!pa)
    
44.                         {
    
45.                                 DeleteKey(keyname,sizeof(luVOID),Matrix,DelMatrix,1);        //将矩阵对象放回缓冲区
    $ A, V" `% Y* C5 {
46.                                 return NULL;
    6 J" r/ L: t/ b$ k
47.                         }
    * f/ P( n- c6 u) s/ H! A
48.                         delete[] pMatrix->Array;
    
49.                         pMatrix->Array=pa;
    0 f* j+ _1 p5 x' H; e7 S
50.                 }
    % H6 U$ T8 d; D# s1 _9 B, n3 t
51.         }
    3 S( r6 U, R: S* G* U
52.         else
    % q# O% c+ c" f. E5 [' L
53.         {
    
54.                 pMatrix=new myMatrix;        //创建矩阵对象
    ! {8 f( s: D# D
55.                 if(!pMatrix) return NULL;
    
56.                 pMatrix->Array=new double[k];
    # o3 ?+ }7 n5 d. p2 d; S
57.                 if(!pMatrix->Array)
    
58.                 {
    , i0 `3 x; n. I& _
59.                         delete pMatrix;
    ; v4 J: R' Y8 @3 C: K
60.                         return NULL;
    / [( z9 C) f+ d9 D1 Z* N
61.                 }
    , {, i: B  }* F
62.                 if(InsertKey((char *)&pMatrix,-1,Matrix,pMatrix,DelMatrix,NULL,0,NowKey))        //将矩阵对象注册到Lu
    # h0 v# ]8 Y& i" o8 _
63.                 {
    : a) V5 g3 X3 p6 f
64.                         delete pMatrix;
    ! T  p# G, ?% X  T
65.                         return NULL;
    
66.                 }
    
67.         }
    
68.         pMatrix->ArrayLen=k; pMatrix->Dim[0]=m; pMatrix->Dim[1]=n;
    
69.         return pMatrix;
    % w2 B% ]! q9 v
70. }
    
71. LuData _stdcall OpMatrix(luINT mm,LuData *xx,void *hFor,int theOperator)        //运算符重载函数，用于LockKey函数
    4 d' f+ }& |9 S0 @7 `/ ~
72. {
    
73.         LuData a;
    
74.         myMatrix *pMatrix1,*pMatrix2,*pMatrix3;
    $ r+ ]) \0 a/ N7 e
75.         luVOID i,j,k,m,n,u,v;
    0 t" n- [9 x2 x/ N( N/ l( c$ P
76.         double *pa,*pb,*pc;
    : z) B0 \, O9 ~+ E& V
77.         luMessage pMessage;
    
78.         wchar_t wchNum[32];
    ! ~' o% `+ W5 V$ n( g" E
79.         char chNum[32];
    : {6 x1 R! C) u% z' i, j- n* }" u3 o
80.         a.BType=luStaData_nil; a.VType=luStaData_nil; a.x=0;
    5 o2 Z: ~* b0 B
81.         switch(theOperator)
    
82.         {
    
83.         case 2:        //重载运算符*
    
84.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
    
85.                 pMatrix2=(myMatrix *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),Matrix);
    ! |, a4 V, e0 Q8 S& F9 e
86.                 if(!pMatrix1 || !pMatrix2) break;        //对象句柄无效，不是矩阵
    
87.                 if(pMatrix1->Dim[1]!=pMatrix2->Dim[0]) break;        //维数不匹配
    ( [) Y: i8 r0 E& v: m( S2 [
88.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix2->Dim[1]);        //生成新矩阵
    
89.                 if(!pMatrix3) break;
    : @3 c# `  m7 U  ?! O
90.                 pa=pMatrix1->Array; pb=pMatrix2->Array; pc=pMatrix3->Array;
    
91.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=pMatrix2->Dim[1];
    / M. H& v; u; K. s; F' z# C$ z$ u4 d
92.                 for(i=0; i<m; i++)        //矩阵乘
    
93.                 {
    1 i8 q7 G; i2 b1 f
94.                         for(j=0; j<k; j++)
    
95.                         {
    
96.                                 u=i*k+j; pc[u]=0.0;
    9 i3 E& q9 ~* m( t+ i) {
97.                                 for (v=0; v<n; v++)
    * |+ G" Z( `* g% v0 f5 g
98.                                 {
    
99.                                         pc[u]=pc[u]+pa[i*n+v]*pb[v*k+j];
    
100.                                 }
     / p  W, l8 D* e4 l) ]
101.                         }
     
102.                 }
     
103.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     * D, y2 ~* R# o+ Q" g+ G+ o
104.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     2 H5 x3 e: J  R; V) u
105.                 break;
     ' J& o0 z- o$ v0 V( K% B1 y
106.         case 25:        //重载运算符.*
     , ^/ ]( b: w+ A/ e
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;        //对象句柄无效，不是矩阵
     - Z1 g& K/ K- k) a0 Q% w
110.                 if(pMatrix1->Dim[0]!=pMatrix2->Dim[0] || pMatrix1->Dim[1]!=pMatrix2->Dim[1]) break;        //维数不相同
     % L) `2 E$ s" C" W9 V
111.                 pMatrix3=NewMatrix(pMatrix1->Dim[0],pMatrix1->Dim[1]);        //生成新矩阵
     
112.                 if(!pMatrix3) break;
     * k/ ^9 [) e  k: @  M
113.                 for(i=0;i<pMatrix1->ArrayLen;i++) pMatrix3->Array[i]=pMatrix1->Array[i]*pMatrix2->Array[i];        //矩阵点乘
     8 ]3 ?( ^4 d5 D  k4 o4 f! c" u9 J
114.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
     
115.                 a.BType=Matrix; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pMatrix3;
     
116.                 break;
     
117.         case 46:        //重载函数new
     
118.                 if(mm<2) break;
     # m& t# s5 \( @
119.                 if((xx+1)->x<1 || (xx+2)->x<1 || (xx+1)->BType!=luStaData_int64 || (xx+2)->BType!=luStaData_int64) break;
     5 b  U% P1 y5 W1 D' ?# X
120.                 pMatrix3=NewMatrix((luVOID)(xx+1)->x,(luVOID)(xx+2)->x);//生成新矩阵
     
121.                 if(!pMatrix3) break;
     
122.                 for(j=0,i=3;i<=mm;i++,j++)        //赋初值
     6 h7 T$ `' A+ L5 z# q
123.                 {
     , ?& l* a7 j  `) u. h
124.                         if(j>=pMatrix3->ArrayLen) break;
     ) L, d+ U, V1 ?, \7 [! S+ M5 X
125.                         if((xx+i)->BType!=luStaData_double) break;        //只接受实数参数
     
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;
     % v; y" ?9 a0 \0 k
130.                 break;
     1 L( M2 H6 ^2 Q6 A) z
131.         case 49:        //重载函数o
     
132.                 pMessage=(luMessage)SearchKey("\0\0\0\0",sizeof(luVOID),luPubKey_User);
     
133.                 if(!pMessage) break;
     9 t  i2 H3 T8 h& q5 c
134.                 pMatrix1=(myMatrix *)SearchKey((char *)&(xx->x),sizeof(luVOID),Matrix);
     
135.                 if(!pMatrix1) break;        //对象句柄无效，不是矩阵
     3 E' f0 r; @+ f: q0 k
136.                 pa=pMatrix1->Array;
     
137.                 m=pMatrix1->Dim[0]; n=pMatrix1->Dim[1]; k=0;
     
138.                 for(i=0; i<m; i++)        //输出矩阵
     $ {; W, F7 {6 `2 Y6 ]
139.                 {
     ; j) G0 h, |& i$ W, x
140.                         pMessage(L"\r\n"); k+=2;
     
141.                         for(j=0; j<n; j++)
     
142.                         {
     - C9 U' X! ]0 i2 S
143.                                 _gcvt_s(chNum,pa[i*n+j],16);
     
144.                                 for(u=0;chNum[u];u++) {wchNum[u]=chNum[u]; k++;}
     
145.                                 wchNum[u]='\0';
     
146.                                 pMessage(wchNum); pMessage(L"  "); k+=2;
     
147.                         }
     
148.                 }
     
149.                 pMessage(L"\r\n"); k+=2;
     . X8 T. U1 g9 @( t. m" N
150.                 a.BType=luStaData_int64; a.VType=luStaData_int64; a.x=k;        //按函数o的要求，返回输出的字符总数
     % J( S5 a! y/ P9 b
151.                 break;
     ) G4 Z9 }  V5 w* l2 K  W
152.         default:
     * z, b0 b. O9 T" z$ x" I; g6 {
153.                 break;
     
154.         }
     
155.         return a;
     / _, c- v9 Q/ M
156. }
     
157. void main(void)
     1 }9 i" r. S; w1 \7 a7 W1 j2 o
158. {
     
159.         void *hFor;                //表达式句柄
     , v& ^/ n8 N( R, @# C
160.         luINT nPara;                //存放表达式的自变量个数
     
161.         LuData *pPara;                //存放输入自变量的数组指针
     - s, A( F% y. f% m( ?8 ~
162.         luINT ErrBegin,ErrEnd;        //表达式编译出错的初始位置和结束位置
     $ @7 c1 ]1 k) \+ m, P
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.]}";//字符串表达式，矩阵乘
     ) ]# [6 Y, b& z& a% ~  u
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.]}";//字符串表达式，矩阵点乘
     - U; Y8 t) o6 C# h) e8 r; A- x
167.         LuData Val;
     + e6 M+ s$ T6 ^7 W2 K* L
168.         if(!InitLu()) return;        //初始化Lu
     
169.         while(LockKey(Matrix,DelMatrix,OpMatrix)){Matrix--;}        //锁定一个键，用于存储矩阵扩展类型
     8 ?3 F8 X9 d) v0 ~4 k4 C6 O
170. 
     0 J( G5 ]* D' y7 ]# u& k
171.         Val.BType=luStaData_int64; Val.VType=luStaData_int64; Val.x=Matrix;        //定义整数常量
     , |& b0 z1 Z' s, I" J
172.         SetConst(L"matrix",&Val);        //设置整数常量
     " `) k' A9 [7 \% n% X' A, h- Y
173.         InsertKey("\0\0\0\0",4,luPubKey_User,LuMessage,NULL,NULL,1,v); //使Lu运行时可输出函数信息
     
174.         wcout.imbue(locale("chs"));        //设置输出的locale为中文
     & h1 `. g3 `% g! I$ e1 c
175.   
     
176.         ErrCode=LuCom(ForStr,0,0,0,hFor,nPara,pPara,ErrBegin,ErrEnd); //编译表达式
     ) Y6 }: G6 O% b
177.         if(ErrCode)
     ; g/ i' M, ?; V4 o1 J' `
178.         {
     
179.                 wcout<<L"表达式有错误！错误代码："<<ErrCode<<endl;
     
180.         }
     
181.         else
     " }" J( V$ Y! |5 D5 t" I4 M$ p+ a1 Z
182.         {
     
183.                 LuCal(hFor,pPara);        //计算表达式的值
     
184.         }
     ) u$ [9 U( _' [. F
185.         LockKey(Matrix,NULL,OpMatrix);//解锁键Matrix，本例中，该函数可以不用
     8 r$ N5 D4 c/ h: `! S: @
186.         FreeLu();                        //释放Lu
     
187. }

复制代码

习题：
- J5 D5 p8 _6 T+ Q4 x! i: ^
    （1）自定义矩阵的加、减、左除、右除、点左除等运算，自编测试字符串代码，重新编译运行程序，观察计算结果。

5 L  v' U" E, s6 z; }6 T# D    （2）小矩阵乘效率测试。编译运行以下Lu字符串代码：

1. main(:a,b,c,d,t,i)=
   
2.     a=new[matrix,2,2: 1.,2.,2.,1.],
   " Z+ R( r7 }' j  ~! T: V. ^+ z, G
3.     b=new[matrix,2,2: 2.,1.,1.,2.],
   ! o" c7 |9 T" {% s
4.     c=new[matrix,2,2: 2/3.,-1/3.,-1/3.,2/3.],
   " V, H, v/ L5 W# o, U' `
5.     t=clock(),
   
6.     d=a*b, i=0, while{i<1000000, d=d*c*b, i++},
   
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.];
   : s5 q; `& q  \! @, l, l4 C
2. b=[2.,1.;1.,2.];
   ' d  k3 F0 y' A7 j
3. c=[2/3.,-1/3.;-1/3.,2/3.];
   
4. tic,
   
5. d=a*b;
   / W6 r" |- g* o0 Y  w3 b
6. for i=1:1000000
   / s0 p3 t5 U& w8 F2 U" c
7.     d=d*c*b;
   : K3 k$ \7 |9 V- t& ^
8. end
   
9. d,
   ( e. J. o; u% O) w8 e7 ?
10. toc

复制代码

结果：

1. d =
   
2.      4     5
   % R5 s9 I! N9 t& O1 N
3.      5     4
   
4. Elapsed time is 2.903034 seconds.

复制代码

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

; a6 ]/ B0 |3 r/ I    由以上可以看出，自定义数据类型和系统内置类型有近乎相同的效率。

---

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

Powered by Discuz! X2.5