数学建模社区-数学中国

标题: Lu基于系统内置对象创建扩展数据类型，小矩阵乘效率测试 [打印本页]

---

作者: forcal    时间: 2011-10-23 09:37
标题: Lu基于系统内置对象创建扩展数据类型，小矩阵乘效率测试
    本例中，我们将基于系统内置实数数组创建矩阵（matrix）类型，即：基本类型为luDynData_realarray（标识实数数组），扩展类型为matrix（标识矩阵）。为了简单，我们仅处理二维实数数组即矩阵类型。
) a! W1 [0 C  ?' S0 [4 [9 V) j
    基本要点：

% g$ U; @. c  @4 s0 ^% p+ ]    （1）为扩展类型matrix编写运算符重载函数OpMatrix。

    （2）用函数LockKey将重载函数OpMatrix注册到Lu，锁定的键的类型即为matrix，要注册为常量，以便于使用。

$ h# C* ^! W# ]8 o( P& F% {4 b    （3）为扩展类型matrix编写其他操作函数（本例未提供）。

9 ], x0 `! {  w" b( I' N3 z+ B7 t    （4）用函数LockKey解锁键matrix（本例中，程序退出时会自动解锁，故可以不用）。

1. #include <windows.h>
   - ?7 G& i- e) P- j/ J0 `3 Y9 I7 o7 m
2. #include <iostream>
   
3. #include <math.h>
   
4. #include "lu32.h"
   " v9 i, A  E. T+ Y; p* \. ^
5. 
   3 P% m1 f& }7 ^
6. #pragma comment( lib, "lu32.lib" )
   
7. 
   3 c) \1 l2 }+ q5 P9 K$ h. @
8. using namespace std;
   , [+ b  B5 t% ~+ p/ D# U
9. 
   
10. luKEY Matrix=-1000;        //标识矩阵类型，最终的Matrix由LockKey决定
    / r7 ]0 Z0 d9 A
11. 
    
12. void _stdcall LuMessage(wchar_t *pch)//输出动态库信息，该函数注册到Lu，由Lu二级函数调用
    
13. {
    7 f1 l1 ^: U/ Y6 C+ [
14.         wcout<<pch;
    0 L* m( h& [9 }2 G% t: B
15. }
    
16. void _stdcall DelMatrix(void *me)        //用于LockKey函数，因为是基于系统内置实数数组创建矩阵，故该函数什么也不做
    . G1 F& q  k- y, @0 v) c
17. {
    
18. }
    
19. LuData _stdcall OpMatrix(luINT mm,LuData *xx,void *hFor,int theOperator)        //运算符重载函数，用于LockKey函数
    : W* M% t4 C  A8 J" X0 f
20. {
    ( ~/ i4 y; d! T  |/ j1 ]
21.         LuData a;
    
22.         luRealArray *pRealArray1,*pRealArray2,*pRealArray3;
    
23.         luVOID i,j,k,m,n,u,v;
    $ o5 @# @: T3 Y1 t  w; @
24.         double *pa,*pb,*pc;
    
25.         luMessage pMessage;
    
26.         wchar_t wchNum[32];
    * [! ^) D% m, o6 a1 U
27.         char chNum[32];
    
28.         a.BType=luStaData_nil; a.VType=luStaData_nil; a.x=0;
    
29.         switch(theOperator)
    
30.         {
    5 }3 E9 K" B' F" ]
31.         case 2:        //重载运算符*
    + F+ |+ q& x' N7 W
32.                 pRealArray1=(luRealArray *)SearchKey((char *)&(xx->x),sizeof(luVOID),luDynData_realarray);
    4 ?  ]6 J  x7 D, o* T
33.                 pRealArray2=(luRealArray *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),luDynData_realarray);
    ( y! |$ g( K( i# m: H
34.                 if(!pRealArray1 || !pRealArray2) break;        //对象句柄无效，不是实数数组
    ) E3 L& c; C( y; `8 q" t
35.                 if(pRealArray1->DimLen!=2 || pRealArray2->DimLen!=2) break;        //不是二维实数数组（矩阵）
    ! R2 N  p1 ]; }+ G7 Y6 I  R
36.                 if(pRealArray1->Dim[1]!=pRealArray2->Dim[0]) break;        //维数不匹配
    # z  }$ {& l! l; V2 z: j
37.                 pRealArray3=(luRealArray *)NewSysObj(luDynData_realarray,pRealArray1->Dim[0]*pRealArray2->Dim[1],2);        //创建矩阵对象
    & ^$ ?6 D% z+ i+ W; h6 h8 a8 p
38.                 if(!pRealArray3) break;
    ' f: [$ g" Z  @4 F0 c2 e
39.                 pRealArray3->Dim[0]=pRealArray1->Dim[0]; pRealArray3->Dim[1]=pRealArray2->Dim[1];        //设置矩阵维数大小
    
40.                 pa=pRealArray1->Array; pb=pRealArray2->Array; pc=pRealArray3->Array;
    2 Z' j: r1 t8 e+ o9 g+ t) n
41.                 m=pRealArray1->Dim[0]; n=pRealArray1->Dim[1]; k=pRealArray2->Dim[1];
    
42.                 for(i=0; i<m; i++)        //矩阵乘
      n! o3 l1 B3 g! r
43.                 {
    * I% \0 _0 U# F* b: R6 R; x
44.                         for(j=0; j<k; j++)
    
45.                         {
    
46.                                 u=i*k+j; pc[u]=0.0;
    
47.                                 for (v=0; v<n; v++)
    + H1 d; W+ C- e/ A' s7 C
48.                                 {
    
49.                                         pc[u]=pc[u]+pa[i*n+v]*pb[v*k+j];
    0 U2 I" q6 D6 }" I
50.                                 }
    
51.                         }
    
52.                 }
    ' a! U7 }5 k/ a" q, F( L
53.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
    ; E! ~. _2 ?3 U7 E  c" ^+ [
54.                 a.BType=luDynData_realarray; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pRealArray3;
    % h+ N$ S/ T0 Q6 N. f/ I
55.                 break;
    
56.         case 25:        //重载运算符.*
    6 W( P( U5 o( b
57.                 pRealArray1=(luRealArray *)SearchKey((char *)&(xx->x),sizeof(luVOID),luDynData_realarray);
    
58.                 pRealArray2=(luRealArray *)SearchKey((char *)&((xx+1)->x),sizeof(luVOID),luDynData_realarray);
    ( h( M3 S' T1 ^) ]9 J1 _9 N% ]
59.                 if(!pRealArray1 || !pRealArray2) break;        //对象句柄无效，不是实数数组
    4 {/ A3 P/ `; E2 [& z
60.                 if(pRealArray1->DimLen!=2 || pRealArray2->DimLen!=2) break;        //不是二维实数数组（矩阵）
    
61.                 if(pRealArray1->Dim[0]!=pRealArray2->Dim[0] || pRealArray1->Dim[1]!=pRealArray2->Dim[1]) break;        //维数不相同
    
62.                 pRealArray3=(luRealArray *)NewSysObj(luDynData_realarray,pRealArray1->ArrayLen,2);        //创建矩阵对象
    
63.                 if(!pRealArray3) break;
    
64.                 pRealArray3->Dim[0]=pRealArray1->Dim[0]; pRealArray3->Dim[1]=pRealArray1->Dim[1];        //设置矩阵维数大小
    
65.                 for(i=0;i<pRealArray1->ArrayLen;i++) pRealArray3->Array[i]=pRealArray1->Array[i]*pRealArray2->Array[i];//矩阵点乘
    
66.                 FunReObj(hFor);        //告诉Lu，返回一个动态对象
    
67.                 a.BType=luDynData_realarray; a.VType=Matrix; a.x=0; *(luVOID *)&(a.x)=(luVOID)pRealArray3;
    , w) L& b6 ]9 M7 P
68.                 break;
    
69.         case 46:        //重载函数new
    
70.                 a=ExeOperator(mm,xx,hFor,theOperator,luDynData_realarray);        //直接调用基本类型luDynData_realarray的new函数
    
71.                 if(a.VType==luDynData_realarray) a.VType=Matrix;        //设置扩展类型为自定义的Matrix类型
    
72.                 break;
    , d$ x5 `( u) `& @% A
73.         case 49:        //重载函数o
    
74.                 pMessage=(luMessage)SearchKey("\0\0\0\0",sizeof(luVOID),luPubKey_User);
    : j0 G# r/ ~+ u6 H$ E
75.                 if(!pMessage) break;
    6 K& F* m; A6 g3 }9 x0 v" n
76.                 pRealArray1=(luRealArray *)SearchKey((char *)&(xx->x),sizeof(luVOID),luDynData_realarray);
    
77.                 if(!pRealArray1) break;        //对象句柄无效，不是实数数组
    / S) {5 X  u0 f6 `, B
78.                 if(pRealArray1->DimLen!=2) break;        //不是二维实数数组（矩阵）
    
79.                 pa=pRealArray1->Array;
    
80.                 m=pRealArray1->Dim[0]; n=pRealArray1->Dim[1]; k=0;
    ! e) g4 e" q0 ]; c' C! B  {
81.                 for(i=0; i<m; i++)        //输出矩阵
      V5 M# `/ X; W# J4 |
82.                 {
    5 Q8 v% }4 V7 S7 X- X* Z
83.                         pMessage(L"\r\n"); k+=2;
    
84.                         for(j=0; j<n; j++)
    
85.                         {
    ! s( C( d% f/ O$ e
86.                                 _gcvt_s(chNum,pa[i*n+j],16);
    
87.                                 for(u=0;chNum[u];u++) {wchNum[u]=chNum[u]; k++;}
    ' A0 k2 v9 l) [" C! ~4 c, H& [
88.                                 wchNum[u]='\0';
      _, y% M* W9 K: j4 S/ O0 Z6 s
89.                                 pMessage(wchNum); pMessage(L"  "); k+=2;
    
90.                         }
    
91.                 }
    
92.                 pMessage(L"\r\n"); k+=2;
    ; l- i3 J* t/ t# x$ U7 G1 b
93.                 a.BType=luStaData_int64; a.VType=luStaData_int64; a.x=k;        //按函数o的要求，返回输出的字符总数
    ! o! G! \; E) }# N7 I
94.                 break;
    % e% b, O# _1 T3 h/ R' S
95.         default:
    $ a% g' k  U, c" S) `. A( Y2 P
96.                 break;
    / {; ]1 p; v! r& e0 x
97.         }
    ! X5 s' F$ `& A& x0 j! g
98.         return a;
    
99. }
    
100. void main(void)
     
101. {
     
102.         void *hFor;                //表达式句柄
     
103.         luINT nPara;                //存放表达式的自变量个数
     
104.         LuData *pPara;                //存放输入自变量的数组指针
     : d( T7 D" U+ L" x# w
105.         luINT ErrBegin,ErrEnd;        //表达式编译出错的初始位置和结束位置
     
106.         int ErrCode;                //错误代码
     / r, {0 b% N+ n$ w9 d
107.         void *v;
     , P- H/ g4 L" t) T2 A9 X
108.         wchar_t ForStr[]=L"o{new[matrix,2,3,data: 0.,1.,2.;3.,4.,5.]*new[matrix,3,2,data: 1.,2.;3.,4.;5.,6.]}";//字符串表达式，矩阵乘
     
109.         //wchar_t ForStr[]=L"o{new[matrix,2,3,data: 0.,1.,2.;3.,4.,5.].*new[matrix,2,3,data: 1.,2.,3.;4.,5.,6.]}";//字符串表达式，矩阵点乘
     
110.         LuData Val;
     ' @2 c) c8 e) l' K
111.         if(!InitLu()) return;        //初始化Lu
     + Y- _$ B. w8 p' o4 c- U/ y* n
112.         while(LockKey(Matrix,DelMatrix,OpMatrix)){Matrix++;}        //锁定一个键，用于存储矩阵扩展类型
     6 e7 l0 g0 u6 z1 s( l% q. t- w: N4 J. |
113. 
     ! k  ~4 U: U6 T% A" ?) Z/ ?8 w" V
114.         Val.BType=luStaData_int64; Val.VType=luStaData_int64; Val.x=Matrix;        //定义整数常量
     6 s! Y1 _! V! _. G* k0 y
115.         SetConst(L"matrix",&Val);        //设置整数常量
     
116.         InsertKey("\0\0\0\0",4,luPubKey_User,LuMessage,NULL,NULL,1,v); //使Lu运行时可输出函数信息
     1 b) Q7 X- j2 }$ v
117.         wcout.imbue(locale("chs"));        //设置输出的locale为中文
     
118.   
     
119.         ErrCode=LuCom(ForStr,0,0,0,hFor,nPara,pPara,ErrBegin,ErrEnd); //编译表达式
     $ n, o/ l: c/ C* Z* U+ l% ?) P
120.         if(ErrCode)
     6 j- S: w6 v  Y/ N- f0 h
121.         {
     ' Q" M% e8 T2 _
122.                 wcout<<L"表达式有错误！错误代码："<<ErrCode<<endl;
     ( X! q2 q7 ~0 A% P
123.         }
     
124.         else
     
125.         {
     
126.                 LuCal(hFor,pPara);        //计算表达式的值
     - p. b) M2 S+ I
127.         }
     # F5 z* k5 r8 j+ R4 X0 M( k7 r6 X
128.         LockKey(Matrix,NULL,OpMatrix);//解锁键Matrix，本例中，该函数可以不用
     
129.         FreeLu();                        //释放Lu
     
130. }

复制代码

习题：
5 ]. L! }5 @' u% ?( w( U0 s3 z
$ h1 S# {. `. m6 X) k    （1）自定义矩阵的加、减、左除、右除、点左除等运算，自编测试字符串代码，重新编译运行程序，观察计算结果。
- M5 ?( g9 m) J% M
% z1 X- c9 O' p/ ~+ Y  m    （2）小矩阵乘效率测试。编译运行以下Lu字符串代码：

1. main(:a,b,c,d,t,i)=
   
2.     a=new[matrix,2,2,data:1.,2.,2.,1.],
   , g3 s# V" U! g2 B4 a- \
3.     b=new[matrix,2,2,data:2.,1.,1.,2.],
   : X. C& j. j7 n+ ^
4.     c=new[matrix,2,2,data: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++},
   
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,data:1.1,2.,2.,1.], b=new[matrix,2,2,data:2.,1.,1.,2.], c=new[matrix,2,2,data: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.
   7 W! D; Z& Y1 Y" J5 g% g' K, X
2. 5. 4.
   / e* h2 F$ O9 [+ ?9 v/ n3 U
3. time=0.875 seconds.
     P. o0 \& c2 F; d
4. 请按任意键继续. . .

复制代码

Matlab 2009a 代码：

1. a=[1.,2.;2.,1.];
   ; A6 Y' {4 B0 Y) R. D
2. b=[2.,1.;1.,2.];
   
3. c=[2/3.,-1/3.;-1/3.,2/3.];
   2 p% E6 P& D1 e5 R- n) U
4. tic,
   
5. d=a*b;
   
6. for i=1:1000000
   
7.     d=d*c*b;
   
8. end
   
9. d,
   2 l, {3 F1 j9 O! R, R: @7 j) V
10. toc

复制代码

结果：

1. d =
   
2.      4     5
   : h. N# i4 R% }/ v
3.      5     4
   . q- V, H0 L" E; n4 {" G
4. Elapsed time is 2.903034 seconds.

复制代码

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

---

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

Powered by Discuz! X2.5