-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathmex_lf_c.m
147 lines (134 loc) · 5.09 KB
/
mex_lf_c.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
function mexfuncname=mex_lf_c(pathname,filename,Yred)
%Creates C_code from Maple equation file, filename, that compiles
%as a MEX-function. mexfilename is the name of the C_code file.
%The inputs: filename and pathname can be obtained using uigetfile
%
clear mex;
in_filename='model.src';
mexfilename=filename;
mexfuncname=mexfilename(1:length(mexfilename)-2);
%---------------------------------------
t=clock;
procread('model.src');
maple('model(0):')
'load maple model'
etime(clock,t)
procread('mclassic.src');
t=clock;
dummy_var=0;
[zrow,zcol,Sred]=find(Yred);
F=maple('mclassic',zrow,zcol);
'classical model'
etime(clock,t)
procread('mclasord.src');
procread('mc_subop.src');
ExpLimit=5;
t=clock;
dimensions=double(sym(maple('mclasord',ExpLimit)));
'reordered model'
etime(clock,t)
t=clock;
lengths=double(sym(maple('mc_subop',ExpLimit,dimensions(1),dimensions(2),dimensions(3))));
'suboptimized model'
etime(clock,t)
fid=fopen([pathname,mexfilename],'wt');
%------------------------- write header of C program -------------
fprintf(fid,'/*\n');
fprintf(fid,' The calling syntax is:\n');
fprintf(fid,[' [f,J] = ',mexfuncname,'(data,x,param)\n']);
fprintf(fid,'*/\n\n');
fprintf(fid,'#include <math.h>\n');
fprintf(fid,'#include "mex.h"\n\n');
fprintf(fid,'/* Input Arguments */\n\n');
fprintf(fid,'#define Data_IN prhs[0]\n');
fprintf(fid,'#define X_IN prhs[1]\n');
fprintf(fid,'#define Param_IN prhs[2]\n');
fprintf(fid,'/* Output Arguments */\n\n');
fprintf(fid,'#define F_OUT plhs[0]\n');
fprintf(fid,'#define J_OUT plhs[1]\n\n');
fprintf(fid,'static\n');
fprintf(fid,'#ifdef __STDC__\n');
fprintf(fid,['void ',mexfuncname,'(\n']);
fprintf(fid,' double f[],\n');
fprintf(fid,' double J[],\n');
fprintf(fid,' double data[],\n');
fprintf(fid,' double x[],\n');
fprintf(fid,' double param[]\n');
fprintf(fid,' )\n');
fprintf(fid,'#else\n');
fprintf(fid,[mexfuncname,'(f,J,data,x,param)\n']);
fprintf(fid,'double f[],J[];\n');
fprintf(fid,'double data[],x[],param[];\n');
fprintf(fid,'#endif\n\n');
fprintf(fid,'/* Computational Routine */\n');
fprintf(fid,'{\n');
%-------------- write variable declaration and optimized C code ------
remainder=lengths(2);
while remainder>0
tv=maple(['seq(templist[ii+',num2str(lengths(2)-remainder,8),'],ii=1..min(10,',num2str(remainder,8),'));']);
fprintf(fid,' double %s;\n',tv);
remainder=remainder-10;
end;
fprintf(fid,'\n');
fprintf(fid,' double %s;\n','s1, s2, s3, s4, s5, s6, s7, s8, s9, s10');
fprintf(fid,'\n');
fclose(fid);
%------- write C code ---------------
maple('gc:');
maple('readlib(C):');
'C lib installed'
maple('interface(screenwidth=8000)');
procread('cwrite.src');
procread('cwriteln.src');
procread('c_conv.src');
t=clock;
fid=fopen([pathname,mexfilename],'a+');
NumExp=10;
remainder=lengths(1)
while remainder>0
nd=fprintf(fid,'%s \n', maplemex(['c_conv(',num2str(min(NumExp,remainder),8),',',num2str(lengths(1)-remainder,8),'):']));
remainder=remainder-NumExp;
end;
fclose(fid);
'write C code'
etime(clock,t)
%------ write the end of computational routine and the begining of gateway routine ---------------------
fid=fopen([pathname,mexfilename],'a+');
fprintf(fid,'\n');
fprintf(fid,' return;\n');
fprintf(fid,'}\n\n');
fprintf(fid,'/* Gateway Routine */\n\n');
fprintf(fid,'#ifdef __STDC__\n');
fprintf(fid,'void mexFunction(\n');
fprintf(fid,' int nlhs,\n');
fprintf(fid,' Matrix *plhs[],\n');
fprintf(fid,' int nrhs,\n');
fprintf(fid,' Matrix *prhs[]\n');
fprintf(fid,' )\n');
fprintf(fid,'#else\n');
fprintf(fid,'mexFunction(nlhs, plhs, nrhs, prhs)\n');
fprintf(fid,'int nlhs, nrhs;\n');
fprintf(fid,'Matrix *plhs[], *prhs[];\n');
fprintf(fid,'#endif\n\n');
fprintf(fid,'{\n');
fprintf(fid,' double *f;\n');
fprintf(fid,' double *J;\n');
fprintf(fid,' double *data;\n');
fprintf(fid,' double *x;\n');
fprintf(fid,' double *param;\n');
fprintf(fid,' /* Create a matrix for the return argument */\n\n');
%------------------ write the output matrix dimension ------------
fprintf(fid,' F_OUT = mxCreateFull(%4i, 1, REAL);\n\n',dimensions(1));
fprintf(fid,' J_OUT = mxCreateFull(%4i,%4i, REAL);\n\n',dimensions(2),dimensions(3));
%------------------ write the end of gateway routine ----------------
fprintf(fid,' /* Assign pointers to the various parameters */\n\n');
fprintf(fid,' f = mxGetPr(F_OUT);\n\n');
fprintf(fid,' J = mxGetPr(J_OUT);\n\n');
fprintf(fid,' data = mxGetPr(Data_IN);\n');
fprintf(fid,' x = mxGetPr(X_IN);\n');
fprintf(fid,' param = mxGetPr(Param_IN);\n');
fprintf(fid,' /* Do the actual computations in a subroutine */\n\n');
fprintf(fid,[' ',mexfuncname,'(f,J,data,x,param);\n']);
fprintf(fid,' return;\n');
fprintf(fid,'}\n');
fclose(fid);