/[escript]/trunk/ripley/src/Assemble_PDE_System_2D.c
ViewVC logotype

Contents of /trunk/ripley/src/Assemble_PDE_System_2D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3641 - (show annotations)
Thu Oct 27 02:16:12 2011 UTC (7 years, 11 months ago) by gross
File MIME type: text/plain
File size: 98644 byte(s)
more work on the assemblage
1 /*******************************************************
2 *
3 * Copyright (c) 2003-2010 by University of Queensland
4 * Earth Systems Science Computational Center (ESSCC)
5 * http://www.uq.edu.au/esscc
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 *******************************************************/
12 /**************************************************************/
13 /* */
14 /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
15 /* */
16 /* -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m and -(X_{k,i})_i + Y_k */
17 /* */
18 /* u is a scalar in a 2D domain. The shape functions for test and solution must be identical */
19 /* */
20 /* Shape of the coefficients: */
21 /* A = p.numEqu x DIM x p.numComp x DIM */
22 /* B = p.numEqu x p.numComp x DIM */
23 /* C = p.numEqu x DIM x p.numComp */
24 /* D = p.numEqu x p.numComp */
25 /* X = p.numEqu x DIM */
26 /* Y = p.numEqu */
27 /* */
28 /**************************************************************/
29 #include "Assemble.h"
30 #include "Util.h"
31 #ifdef _OPENMP
32 #include <omp.h>
33 #endif
34 /**************************************************************/
35 void Finley_Assemble_PDE_System_2D(Finley_Assemble_Parameters p,
36 Finley_ElementFile* elements,
37 Paso_SystemMatrix* Mat, escriptDataC* F,
38 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y) {
39 #define DIM 2
40 index_t color;
41 dim_t e, isub;
42 __const double *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p, *A_q, *B_q, *C_q, *D_q, *X_q, *Y_q;
43 double *EM_S, *EM_F, *Vol, *DSDX;
44 index_t *row_index;
45 register dim_t q, s,r,k,m;
46 register double rtmp, rtmp0, rtmp1, rtmp2, rtmp00, rtmp01, rtmp02, rtmp10, rtmp11, rtmp12, rtmp20, rtmp21, rtmp22;
47 bool_t add_EM_F, add_EM_S;
48 bool_t extendedA=isExpanded(A);
49 bool_t extendedB=isExpanded(B);
50 bool_t extendedC=isExpanded(C);
51 bool_t extendedD=isExpanded(D);
52 bool_t extendedX=isExpanded(X);
53 bool_t extendedY=isExpanded(Y);
54 double *F_p=(requireWrite(F), getSampleDataRW(F,0)); /* use comma, to get around the mixed code and declarations thing */
55 dim_t len_EM_S=p.row_numShapesTotal*p.col_numShapesTotal*p.numEqu*p.numComp;
56 dim_t len_EM_F=p.row_numShapesTotal*p.numEqu;
57 {
58 /* GENERATOR SNIP_PRE TOP */
59 const double w56 = -0.048598662630909980773*h1;
60 const double w13 = -0.049382716049382716049*h1/h0;
61 const double w67 = -0.0034784464623227873914*h0;
62 const double w23 = -0.0077160493827160493827*h0/h1;
63 const double w60 = -0.0061728395061728395062*h0;
64 const double w10 = -0.030864197530864197531*h1/h0;
65 const double w74 = 0.0034784464623227873914*h1;
66 const double w95 = -0.083333333333333333333*h0;
67 const double w127 = 0.10954300427416564056*h1;
68 const double w65 = -0.00086961161558069684786*h0;
69 const double w98 = 0.083333333333333333333*h0;
70 const double w68 = -0.00086961161558069684786*h1;
71 const double w128 = 0.061728395061728395062*h1;
72 const double w130 = 0.0086961161558069684786*h1;
73 const double w75 = 0.00086961161558069684786*h0;
74 const double w91 = 0.048598662630909980773*h0;
75 const double w121 = -0.098765432098765432099*h1;
76 const double w102 = 0.0060748328288637475966*h0*h1;
77 const double w27 = -0.049382716049382716049*h0/h1;
78 const double w29 = 0.060748328288637475966*h0/h1;
79 const double w28 = 0.060748328288637475966*h1/h0;
80 const double w148 = 0.25*h0*h1;
81 const double w107 = 0.00039202670923636763834*h0*h1;
82 const double w86 = -0.048598662630909980773*h0;
83 const double w2 = -0.060748328288637475966;
84 const double w45 = 0.25000000000000000000;
85 const double w55 = -0.0068464377671353525349*h0;
86 const double w12 = 0.012345679012345679012*h0/h1;
87 const double w25 = -0.030864197530864197531*h0/h1;
88 const double w8 = 0.060748328288637475966;
89 const double w62 = -0.024691358024691358025*h0;
90 const double w125 = -0.013913785849291149566*h1;
91 const double w104 = 0.0030864197530864197531*h0*h1;
92 const double w109 = 0.047827057692638375834*h0*h1;
93 const double w48 = -0.16666666666666666667*h1/h0;
94 const double w40 = 0.012345679012345679012*h1/h0;
95 const double w119 = -0.061728395061728395062*h1;
96 const double w32 = 0.030864197530864197531*h1/h0;
97 const double w81 = 0.053901890521502123431*h1;
98 const double w19 = 0.00098006677309091909584;
99 const double w39 = -0.00098006677309091909584*h0/h1;
100 const double w77 = 0.0034784464623227873914*h0;
101 const double w14 = 0.049382716049382716049;
102 const double w138 = -0.5*h1;
103 const double w110 = 0.0000124484802011303384*h0*h1;
104 const double w1 = 0.0077160493827160493827;
105 const double w42 = -0.0015681068369454705533*h0/h1;
106 const double w99 = -0.16666666666666666667*h0;
107 const double w115 = -0.068464377671353525349*h0;
108 const double w58 = -0.053901890521502123431*h1;
109 const double w118 = -0.0086961161558069684786*h0;
110 const double w103 = 0.024299331315454990386*h0*h1;
111 const double w140 = 0.5*h1;
112 const double w22 = -0.0077160493827160493827*h1/h0;
113 const double w97 = 0.16666666666666666667*h1;
114 const double w146 = 0.0069568929246455747829*h0*h1;
115 const double w87 = -0.00078405341847273527667*h0;
116 const double w51 = 0.33333333333333333333*h0/h1;
117 const double w61 = -0.024691358024691358025*h1;
118 const double w83 = 0.048598662630909980773*h1;
119 const double w124 = -0.0086961161558069684786*h1;
120 const double w117 = -0.061728395061728395062*h0;
121 const double w20 = -0.0015681068369454705533*h1/h0;
122 const double w84 = 0.027385751068541410139*h0;
123 const double w96 = 0.083333333333333333333*h1;
124 const double w137 = 0.013913785849291149566*h0;
125 const double w11 = 0.0069568929246455747828;
126 const double w135 = 0.10954300427416564056*h0;
127 const double w26 = -0.012345679012345679012*h0/h1;
128 const double w116 = -0.10954300427416564056*h1;
129 const double w92 = 0.053901890521502123431*h0;
130 const double w38 = 0.0077160493827160493827*h1/h0;
131 const double w17 = -0.0069568929246455747828;
132 const double w69 = 0.0068464377671353525349*h1;
133 const double w113 = 0.11111111111111111111*h0*h1;
134 const double w141 = 0.5*h0;
135 const double w36 = 0.00098006677309091909584*h1/h0;
136 const double w131 = 0.013913785849291149566*h1;
137 const double w50 = 0.33333333333333333333*h1/h0;
138 const double w16 = 0.049382716049382716049*h0/h1;
139 const double w143 = 0.054771502137082820279*h0*h1;
140 const double w0 = -0.060748328288637475966*h1/h0;
141 const double w144 = 0.0077160493827160493827*h0*h1;
142 const double w85 = -0.053901890521502123431*h0;
143 const double w79 = 0.0061728395061728395062*h0;
144 const double w3 = 0.0077160493827160493827*h0/h1;
145 const double w35 = 0.0015681068369454705533*h0/h1;
146 const double w53 = -0.33333333333333333333*h0/h1;
147 const double w114 = -0.068464377671353525349*h1;
148 const double w49 = -0.16666666666666666667*h0/h1;
149 const double w126 = 0.068464377671353525349*h1;
150 const double w44 = -0.33333333333333333333*h1/h0;
151 const double w4 = -0.097197325261819961546*h1/h0;
152 const double w82 = 0.0068464377671353525349*h0;
153 const double w46 = -0.25000000000000000000;
154 const double w139 = -0.5*h0;
155 const double w94 = -0.16666666666666666667*h1;
156 const double w37 = 0.0015681068369454705533*h1/h0;
157 const double w64 = -0.00011045515751022224798*h1;
158 const double w78 = 0.00011045515751022224798*h1;
159 const double w89 = 0.00011045515751022224798*h0;
160 const double w111 = 0.055555555555555555556*h0*h1;
161 const double w101 = -0.083333333333333333333*h1;
162 const double w136 = 0.098765432098765432099*h0;
163 const double w9 = -0.0077160493827160493827;
164 const double w31 = 0.00098006677309091909584*h0/h1;
165 const double w18 = -0.00098006677309091909584*h1/h0;
166 const double w73 = 0.024691358024691358025*h1;
167 const double w132 = 0.068464377671353525349*h0;
168 const double w122 = -0.098765432098765432099*h0;
169 const double w34 = 0.049382716049382716049*h1/h0;
170 const double w76 = 0.00078405341847273527667*h1;
171 const double w108 = 0.00077160493827160493827*h0*h1;
172 const double w54 = -0.0068464377671353525349*h1;
173 const double w33 = 0.097197325261819961546*h0/h1;
174 const double w47 = 0.16666666666666666667*h0/h1;
175 const double w80 = 0.024691358024691358025*h0;
176 const double w21 = -0.00098006677309091909584;
177 const double w88 = -0.00011045515751022224798*h0;
178 const double w129 = 0.098765432098765432099*h1;
179 const double w63 = -0.027385751068541410139*h1;
180 const double w123 = -0.013913785849291149566*h0;
181 const double w30 = 0.097197325261819961546*h1/h0;
182 const double w90 = 0.00078405341847273527667*h0;
183 const double w105 = 0.012345679012345679012*h0*h1;
184 const double w7 = 0.030864197530864197531*h0/h1;
185 const double w147 = 0.00098006677309091909584*h0*h1;
186 const double w142 = 0.060748328288637475966*h0*h1;
187 const double w120 = -0.10954300427416564056*h0;
188 const double w66 = -0.00078405341847273527667*h1;
189 const double w145 = 0.049382716049382716049*h0*h1;
190 const double w71 = 0.00086961161558069684786*h1;
191 const double w24 = -0.012345679012345679012*h1/h0;
192 const double w43 = -0.097197325261819961546*h0/h1;
193 const double w72 = 0.027385751068541410139*h1;
194 const double w106 = 0.000098006677309091909584*h0*h1;
195 const double w6 = -0.054771502137082820279;
196 const double w41 = -0.060748328288637475966*h0/h1;
197 const double w59 = -0.0034784464623227873914*h1;
198 const double w134 = 0.0086961161558069684786*h0;
199 const double w112 = 0.027777777777777777778*h0*h1;
200 const double w57 = -0.027385751068541410139*h0;
201 const double w52 = 0.16666666666666666667*h1/h0;
202 const double w5 = 0.054771502137082820279;
203 const double w70 = 0.0061728395061728395062*h1;
204 const double w93 = -0.0061728395061728395062*h1;
205 const double w100 = 0.16666666666666666667*h0;
206 const double w133 = 0.061728395061728395062*h0;
207 const double w15 = -0.049382716049382716049;
208 /* GENERATOR SNIP_PRE BOTTOM */
209 #pragma omp parallel private(EM_S, EM_F,k2_0, k0, k1, k2, *A_p, *B_p, *C_p, *D_p, *X_p, *Y_p)
210 {
211 EM_S=THREAD_MEMALLOC(len_EM_S,double);
212 EM_F=THREAD_MEMALLOC(len_EM_F,double);
213 if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
214 for (k2_0 = 0; k2_0 <2; k2_0++) { /* coloring */
215 #pragma omp parallel for private(i2, i1,i0)
216 for (k2 = k2_0; k2< N0; k2=k2+2) {
217 for (k1 = 0; k1< N1; ++k1) {
218 for (k0 = 0; k0< N0; ++k0) {
219 bool_t add_EM_F=FALSE;
220 bool_t add_EM_S=FALSE;
221 index_t e = k0 + M0 * k1 + M0*M1 * k2;
222 A_p=getSampleDataRO(A,e);
223 B_p=getSampleDataRO(B,e);
224 C_p=getSampleDataRO(C,e);
225 D_p=getSampleDataRO(D,e);
226 X_p=getSampleDataRO(X,e);
227 Y_p=getSampleDataRO(Y,e);
228 for (q=0;q<len_EM_S;++q) EM_S[q]=0;
229 for (q=0;q<len_EM_F;++q) EM_F[q]=0;
230 /* GENERATOR SNIP TOP */
231 /**************************************************************/
232 /* process A: */
233 /**************************************************************/
234 if (NULL!=A_p) {
235 add_EM_S=TRUE;
236 if (extendedA) {
237 for (k=0;k<p.numEqu;k++) {
238 for (m=0;m<p.numComp;m++) {
239 const register double A_00_0 = A_p[INDEX5(k,0,m,0,0, p.numEqu,2,p.numComp,2)];
240 const register double A_01_0 = A_p[INDEX5(k,0,m,1,0, p.numEqu,2,p.numComp,2)];
241 const register double A_10_0 = A_p[INDEX5(k,1,m,0,0, p.numEqu,2,p.numComp,2)];
242 const register double A_11_0 = A_p[INDEX5(k,1,m,1,0, p.numEqu,2,p.numComp,2)];
243 const register double A_00_1 = A_p[INDEX5(k,0,m,0,1, p.numEqu,2,p.numComp,2)];
244 const register double A_01_1 = A_p[INDEX5(k,0,m,1,1, p.numEqu,2,p.numComp,2)];
245 const register double A_10_1 = A_p[INDEX5(k,1,m,0,1, p.numEqu,2,p.numComp,2)];
246 const register double A_11_1 = A_p[INDEX5(k,1,m,1,1, p.numEqu,2,p.numComp,2)];
247 const register double A_00_2 = A_p[INDEX5(k,0,m,0,2, p.numEqu,2,p.numComp,2)];
248 const register double A_01_2 = A_p[INDEX5(k,0,m,1,2, p.numEqu,2,p.numComp,2)];
249 const register double A_10_2 = A_p[INDEX5(k,1,m,0,2, p.numEqu,2,p.numComp,2)];
250 const register double A_11_2 = A_p[INDEX5(k,1,m,1,2, p.numEqu,2,p.numComp,2)];
251 const register double A_00_3 = A_p[INDEX5(k,0,m,0,3, p.numEqu,2,p.numComp,2)];
252 const register double A_01_3 = A_p[INDEX5(k,0,m,1,3, p.numEqu,2,p.numComp,2)];
253 const register double A_10_3 = A_p[INDEX5(k,1,m,0,3, p.numEqu,2,p.numComp,2)];
254 const register double A_11_3 = A_p[INDEX5(k,1,m,1,3, p.numEqu,2,p.numComp,2)];
255 const register double A_00_4 = A_p[INDEX5(k,0,m,0,4, p.numEqu,2,p.numComp,2)];
256 const register double A_01_4 = A_p[INDEX5(k,0,m,1,4, p.numEqu,2,p.numComp,2)];
257 const register double A_10_4 = A_p[INDEX5(k,1,m,0,4, p.numEqu,2,p.numComp,2)];
258 const register double A_11_4 = A_p[INDEX5(k,1,m,1,4, p.numEqu,2,p.numComp,2)];
259 const register double A_00_5 = A_p[INDEX5(k,0,m,0,5, p.numEqu,2,p.numComp,2)];
260 const register double A_01_5 = A_p[INDEX5(k,0,m,1,5, p.numEqu,2,p.numComp,2)];
261 const register double A_10_5 = A_p[INDEX5(k,1,m,0,5, p.numEqu,2,p.numComp,2)];
262 const register double A_11_5 = A_p[INDEX5(k,1,m,1,5, p.numEqu,2,p.numComp,2)];
263 const register double A_00_6 = A_p[INDEX5(k,0,m,0,6, p.numEqu,2,p.numComp,2)];
264 const register double A_01_6 = A_p[INDEX5(k,0,m,1,6, p.numEqu,2,p.numComp,2)];
265 const register double A_10_6 = A_p[INDEX5(k,1,m,0,6, p.numEqu,2,p.numComp,2)];
266 const register double A_11_6 = A_p[INDEX5(k,1,m,1,6, p.numEqu,2,p.numComp,2)];
267 const register double A_00_7 = A_p[INDEX5(k,0,m,0,7, p.numEqu,2,p.numComp,2)];
268 const register double A_01_7 = A_p[INDEX5(k,0,m,1,7, p.numEqu,2,p.numComp,2)];
269 const register double A_10_7 = A_p[INDEX5(k,1,m,0,7, p.numEqu,2,p.numComp,2)];
270 const register double A_11_7 = A_p[INDEX5(k,1,m,1,7, p.numEqu,2,p.numComp,2)];
271 const register double A_00_8 = A_p[INDEX5(k,0,m,0,8, p.numEqu,2,p.numComp,2)];
272 const register double A_01_8 = A_p[INDEX5(k,0,m,1,8, p.numEqu,2,p.numComp,2)];
273 const register double A_10_8 = A_p[INDEX5(k,1,m,0,8, p.numEqu,2,p.numComp,2)];
274 const register double A_11_8 = A_p[INDEX5(k,1,m,1,8, p.numEqu,2,p.numComp,2)];
275 const register double tmp34_0 = A_01_5 + A_01_7;
276 const register double tmp14_0 = A_01_5 + A_01_7 + A_10_1 + A_10_3;
277 const register double tmp38_0 = A_01_3 + A_01_7 + A_10_3 + A_10_7;
278 const register double tmp1_0 = A_00_0 + A_00_2;
279 const register double tmp29_0 = A_01_1 + A_01_5 + A_10_3 + A_10_7;
280 const register double tmp8_0 = A_10_2 + A_10_6;
281 const register double tmp27_0 = A_01_6 + A_10_2;
282 const register double tmp19_0 = A_01_1 + A_01_3 + A_10_5 + A_10_7;
283 const register double tmp26_0 = A_01_2 + A_10_6;
284 const register double tmp16_0 = A_00_0 + A_00_2 + A_00_6 + A_00_8;
285 const register double tmp12_0 = A_01_0 + A_10_8;
286 const register double tmp0_0 = A_01_0 + A_01_8;
287 const register double tmp24_0 = A_01_8 + A_10_8;
288 const register double tmp23_0 = A_01_5 + A_01_7 + A_10_5 + A_10_7;
289 const register double tmp21_0 = A_11_0 + A_11_6;
290 const register double tmp32_0 = A_10_1 + A_10_5;
291 const register double tmp17_0 = A_01_8 + A_10_0;
292 const register double tmp35_0 = A_01_2 + A_01_6;
293 const register double tmp20_0 = A_01_0 + A_10_0;
294 const register double tmp11_0 = A_00_6 + A_00_8;
295 const register double tmp7_0 = A_10_1 + A_10_3;
296 const register double tmp37_0 = A_01_6 + A_10_6;
297 const register double tmp18_0 = A_00_1 + A_00_7;
298 const register double tmp10_0 = A_10_5 + A_10_7;
299 const register double tmp9_0 = A_11_3 + A_11_5;
300 const register double tmp5_0 = A_11_1 + A_11_7;
301 const register double tmp30_0 = A_01_0 + A_01_8 + A_10_0 + A_10_8;
302 const register double tmp40_0 = A_01_2 + A_10_2;
303 const register double tmp22_0 = A_11_2 + A_11_8;
304 const register double tmp4_0 = A_01_3 + A_01_7;
305 const register double tmp28_0 = A_01_3 + A_01_7 + A_10_1 + A_10_5;
306 const register double tmp13_0 = A_01_2 + A_01_6 + A_10_2 + A_10_6;
307 const register double tmp33_0 = A_10_3 + A_10_7;
308 const register double tmp6_0 = A_00_3 + A_00_5;
309 const register double tmp15_0 = A_01_4 + A_10_4;
310 const register double tmp3_0 = A_01_1 + A_01_5;
311 const register double tmp25_0 = A_01_1 + A_01_3 + A_10_1 + A_10_3;
312 const register double tmp31_0 = A_10_0 + A_10_8;
313 const register double tmp2_0 = A_11_0 + A_11_2 + A_11_6 + A_11_8;
314 const register double tmp36_0 = A_01_1 + A_01_3;
315 const register double tmp39_0 = A_01_1 + A_01_5 + A_10_1 + A_10_5;
316 const register double tmp29_1 = tmp17_0*w19;
317 const register double tmp100_1 = tmp22_0*w39;
318 const register double tmp115_1 = tmp40_0*w21;
319 const register double tmp89_1 = A_11_3*w42;
320 const register double tmp83_1 = tmp33_0*w11;
321 const register double tmp75_1 = tmp28_0*w6;
322 const register double tmp30_1 = A_11_4*w27;
323 const register double tmp4_1 = tmp3_0*w5;
324 const register double tmp87_1 = A_01_4*w15;
325 const register double tmp93_1 = tmp22_0*w41;
326 const register double tmp81_1 = A_01_8*w2;
327 const register double tmp99_1 = tmp14_0*w5;
328 const register double tmp1_1 = tmp1_0*w0;
329 const register double tmp70_1 = tmp21_0*w31;
330 const register double tmp95_1 = A_10_6*w19;
331 const register double tmp38_1 = tmp3_0*w11;
332 const register double tmp79_1 = tmp21_0*w39;
333 const register double tmp77_1 = tmp30_0*w9;
334 const register double tmp25_1 = tmp14_0*w11;
335 const register double tmp117_1 = tmp26_0*w2;
336 const register double tmp2_1 = tmp2_0*w3;
337 const register double tmp72_1 = tmp15_0*w15;
338 const register double tmp18_1 = A_00_7*w20;
339 const register double tmp44_1 = A_01_2*w19;
340 const register double tmp52_1 = tmp24_0*w19;
341 const register double tmp15_1 = A_01_4*w14;
342 const register double tmp50_1 = tmp23_0*w11;
343 const register double tmp102_1 = tmp33_0*w5;
344 const register double tmp21_1 = tmp11_0*w18;
345 const register double tmp64_1 = tmp25_0*w11;
346 const register double tmp109_1 = tmp21_0*w41;
347 const register double tmp106_1 = A_11_5*w42;
348 const register double tmp35_1 = A_10_8*w2;
349 const register double tmp16_1 = tmp10_0*w17;
350 const register double tmp49_1 = tmp22_0*w31;
351 const register double tmp118_1 = tmp29_0*w6;
352 const register double tmp114_1 = tmp39_0*w17;
353 const register double tmp19_1 = A_10_8*w21;
354 const register double tmp68_1 = A_11_5*w33;
355 const register double tmp73_1 = tmp26_0*w21;
356 const register double tmp46_1 = tmp20_0*w8;
357 const register double tmp36_1 = tmp4_0*w5;
358 const register double tmp12_1 = tmp9_0*w12;
359 const register double tmp111_1 = A_10_2*w19;
360 const register double tmp31_1 = tmp18_0*w24;
361 const register double tmp105_1 = A_10_6*w8;
362 const register double tmp10_1 = A_00_4*w13;
363 const register double tmp82_1 = tmp32_0*w5;
364 const register double tmp51_1 = A_11_5*w35;
365 const register double tmp57_1 = tmp1_0*w28;
366 const register double tmp80_1 = tmp16_0*w38;
367 const register double tmp48_1 = tmp21_0*w29;
368 const register double tmp47_1 = tmp6_0*w32;
369 const register double tmp110_1 = A_01_8*w21;
370 const register double tmp85_1 = tmp35_0*w9;
371 const register double tmp28_1 = tmp16_0*w22;
372 const register double tmp34_1 = tmp11_0*w0;
373 const register double tmp55_1 = A_00_4*w34;
374 const register double tmp123_1 = tmp40_0*w2;
375 const register double tmp91_1 = A_11_5*w43;
376 const register double tmp7_1 = tmp5_0*w7;
377 const register double tmp41_1 = tmp7_0*w17;
378 const register double tmp0_1 = tmp0_0*w1;
379 const register double tmp32_1 = tmp2_0*w23;
380 const register double tmp98_1 = tmp12_0*w19;
381 const register double tmp107_1 = tmp34_0*w17;
382 const register double tmp76_1 = tmp29_0*w17;
383 const register double tmp120_1 = tmp37_0*w21;
384 const register double tmp40_1 = A_01_6*w8;
385 const register double tmp53_1 = A_00_1*w30;
386 const register double tmp116_1 = tmp27_0*w21;
387 const register double tmp101_1 = A_01_0*w2;
388 const register double tmp45_1 = tmp1_0*w18;
389 const register double tmp67_1 = A_00_7*w30;
390 const register double tmp61_1 = A_00_1*w37;
391 const register double tmp17_1 = A_11_4*w16;
392 const register double tmp66_1 = tmp20_0*w19;
393 const register double tmp84_1 = tmp34_0*w6;
394 const register double tmp92_1 = tmp18_0*w40;
395 const register double tmp23_1 = tmp13_0*w1;
396 const register double tmp97_1 = tmp19_0*w11;
397 const register double tmp5_1 = A_00_1*w4;
398 const register double tmp62_1 = tmp24_0*w8;
399 const register double tmp39_1 = tmp10_0*w6;
400 const register double tmp9_1 = tmp7_0*w6;
401 const register double tmp63_1 = tmp1_0*w36;
402 const register double tmp65_1 = A_11_3*w35;
403 const register double tmp69_1 = tmp11_0*w28;
404 const register double tmp96_1 = tmp17_0*w8;
405 const register double tmp108_1 = A_11_3*w43;
406 const register double tmp14_1 = A_10_4*w15;
407 const register double tmp6_1 = tmp4_0*w11;
408 const register double tmp3_1 = A_10_0*w2;
409 const register double tmp22_1 = tmp12_0*w8;
410 const register double tmp26_1 = tmp9_0*w26;
411 const register double tmp8_1 = tmp6_0*w10;
412 const register double tmp113_1 = tmp38_0*w6;
413 const register double tmp33_1 = tmp19_0*w5;
414 const register double tmp60_1 = tmp22_0*w29;
415 const register double tmp27_1 = tmp15_0*w14;
416 const register double tmp54_1 = A_11_3*w33;
417 const register double tmp71_1 = tmp23_0*w5;
418 const register double tmp37_1 = A_00_7*w4;
419 const register double tmp43_1 = A_10_0*w21;
420 const register double tmp59_1 = A_00_7*w37;
421 const register double tmp58_1 = tmp25_0*w5;
422 const register double tmp78_1 = tmp31_0*w1;
423 const register double tmp20_1 = A_01_6*w19;
424 const register double tmp119_1 = tmp28_0*w17;
425 const register double tmp121_1 = tmp39_0*w6;
426 const register double tmp122_1 = tmp38_0*w17;
427 const register double tmp88_1 = A_10_4*w14;
428 const register double tmp56_1 = tmp11_0*w36;
429 const register double tmp103_1 = tmp32_0*w11;
430 const register double tmp104_1 = tmp36_0*w6;
431 const register double tmp94_1 = A_01_0*w21;
432 const register double tmp11_1 = tmp8_0*w9;
433 const register double tmp86_1 = A_10_2*w8;
434 const register double tmp42_1 = A_00_1*w20;
435 const register double tmp74_1 = tmp27_0*w2;
436 const register double tmp24_1 = tmp5_0*w25;
437 const register double tmp112_1 = tmp37_0*w2;
438 const register double tmp90_1 = tmp36_0*w17;
439 const register double tmp13_1 = A_01_2*w8;
440 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp20_1 + tmp21_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
441 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp33_1 + tmp8_1;
442 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp12_1 + tmp14_1 + tmp15_1 + tmp17_1 + tmp2_1 + tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp7_1 + tmp8_1;
443 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp17_1 + tmp23_1 + tmp27_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1 + tmp7_1;
444 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp17_1 + tmp23_1 + tmp27_1 + tmp47_1 + tmp55_1 + tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1 + tmp7_1;
445 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp24_1 + tmp26_1 + tmp28_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp8_1;
446 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp24_1 + tmp30_1 + tmp47_1 + tmp55_1 + tmp78_1 + tmp79_1 + tmp80_1 + tmp81_1 + tmp82_1 + tmp83_1 + tmp84_1 + tmp85_1 + tmp86_1 + tmp87_1 + tmp88_1 + tmp89_1 + tmp90_1 + tmp91_1 + tmp92_1 + tmp93_1 + tmp94_1 + tmp95_1;
447 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp23_1 + tmp24_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp8_1 + tmp96_1 + tmp97_1 + tmp98_1 + tmp99_1;
448 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp100_1 + tmp101_1 + tmp102_1 + tmp103_1 + tmp104_1 + tmp105_1 + tmp106_1 + tmp107_1 + tmp108_1 + tmp109_1 + tmp110_1 + tmp111_1 + tmp24_1 + tmp30_1 + tmp47_1 + tmp55_1 + tmp78_1 + tmp80_1 + tmp85_1 + tmp87_1 + tmp88_1 + tmp92_1;
449 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp100_1 + tmp106_1 + tmp108_1 + tmp109_1 + tmp11_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp19_1 + tmp24_1 + tmp30_1 + tmp36_1 + tmp38_1 + tmp3_1 + tmp40_1 + tmp44_1 + tmp47_1 + tmp55_1 + tmp80_1 + tmp92_1 + tmp9_1;
450 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp11_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp20_1 + tmp24_1 + tmp30_1 + tmp35_1 + tmp39_1 + tmp41_1 + tmp43_1 + tmp47_1 + tmp4_1 + tmp55_1 + tmp6_1 + tmp79_1 + tmp80_1 + tmp89_1 + tmp91_1 + tmp92_1 + tmp93_1;
451 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp102_1 + tmp103_1 + tmp105_1 + tmp10_1 + tmp111_1 + tmp12_1 + tmp17_1 + tmp2_1 + tmp34_1 + tmp37_1 + tmp42_1 + tmp45_1 + tmp78_1 + tmp7_1 + tmp81_1 + tmp84_1 + tmp85_1 + tmp87_1 + tmp88_1 + tmp8_1 + tmp90_1 + tmp94_1;
452 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp112_1 + tmp113_1 + tmp114_1 + tmp115_1 + tmp17_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp51_1 + tmp54_1 + tmp55_1 + tmp61_1 + tmp63_1 + tmp67_1 + tmp69_1 + tmp72_1 + tmp77_1 + tmp7_1;
453 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp101_1 + tmp104_1 + tmp107_1 + tmp10_1 + tmp110_1 + tmp12_1 + tmp17_1 + tmp18_1 + tmp1_1 + tmp21_1 + tmp2_1 + tmp5_1 + tmp78_1 + tmp7_1 + tmp82_1 + tmp83_1 + tmp85_1 + tmp86_1 + tmp87_1 + tmp88_1 + tmp8_1 + tmp95_1;
454 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp116_1 + tmp117_1 + tmp118_1 + tmp119_1 + tmp24_1 + tmp26_1 + tmp28_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp72_1 + tmp77_1 + tmp8_1;
455 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp120_1 + tmp121_1 + tmp122_1 + tmp123_1 + tmp17_1 + tmp47_1 + tmp53_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp59_1 + tmp60_1 + tmp65_1 + tmp68_1 + tmp70_1 + tmp72_1 + tmp77_1 + tmp7_1;
456 }
457 }
458 } else { /* constant data */
459 for (k=0;k<p.numEqu;k++) {
460 for (m=0;m<p.numComp;m++) {
461 const register double A_00 = A_p[INDEX4(k,0,m,0 p.numEqu,2, p.numComp)];
462 const register double A_01 = A_p[INDEX4(k,0,m,1 p.numEqu,2, p.numComp)];
463 const register double A_10 = A_p[INDEX4(k,1,m,0 p.numEqu,2, p.numComp)];
464 const register double A_11 = A_p[INDEX4(k,1,m,1 p.numEqu,2, p.numComp)];
465 const register double tmp0_0 = A_01 + A_10;
466 const register double tmp5_1 = A_00*w48;
467 const register double tmp3_1 = A_01*w45;
468 const register double tmp4_1 = A_11*w49;
469 const register double tmp11_1 = A_00*w52;
470 const register double tmp7_1 = A_00*w50;
471 const register double tmp2_1 = A_10*w46;
472 const register double tmp0_1 = A_00*w44;
473 const register double tmp13_1 = A_10*w45;
474 const register double tmp10_1 = A_01*w46;
475 const register double tmp6_1 = tmp0_0*w45;
476 const register double tmp8_1 = A_11*w51;
477 const register double tmp12_1 = A_11*w53;
478 const register double tmp9_1 = tmp0_0*w46;
479 const register double tmp1_1 = A_11*w47;
480 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
481 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
482 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1;
483 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
484 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp6_1 + tmp7_1 + tmp8_1;
485 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1 + tmp9_1;
486 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
487 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1 + tmp6_1;
488 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1;
489 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp11_1 + tmp12_1 + tmp2_1 + tmp3_1;
490 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp11_1 + tmp12_1 + tmp2_1 + tmp3_1;
491 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp1_1;
492 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
493 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp13_1 + tmp1_1;
494 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1 + tmp9_1;
495 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp7_1 + tmp8_1 + tmp9_1;
496 }
497 }
498 }
499 }
500 /**************************************************************/
501 /* process B: */
502 /**************************************************************/
503 if (NULL!=B_p) {
504 add_EM_S=TRUE;
505 if (extendedB) {
506 for (k=0;k<p.numEqu;k++) {
507 for (m=0;m<p.numComp;m++) {
508 const register double B_0_0 = B_p[INDEX4(k,0,m,0, p.numEqu,2,p.numComp)];
509 const register double B_1_0 = B_p[INDEX4(k,1,m,0, p.numEqu,2,p.numComp)];
510 const register double B_0_1 = B_p[INDEX4(k,0,m,1, p.numEqu,2,p.numComp)];
511 const register double B_1_1 = B_p[INDEX4(k,1,m,1, p.numEqu,2,p.numComp)];
512 const register double B_0_2 = B_p[INDEX4(k,0,m,2, p.numEqu,2,p.numComp)];
513 const register double B_1_2 = B_p[INDEX4(k,1,m,2, p.numEqu,2,p.numComp)];
514 const register double B_0_3 = B_p[INDEX4(k,0,m,3, p.numEqu,2,p.numComp)];
515 const register double B_1_3 = B_p[INDEX4(k,1,m,3, p.numEqu,2,p.numComp)];
516 const register double B_0_4 = B_p[INDEX4(k,0,m,4, p.numEqu,2,p.numComp)];
517 const register double B_1_4 = B_p[INDEX4(k,1,m,4, p.numEqu,2,p.numComp)];
518 const register double B_0_5 = B_p[INDEX4(k,0,m,5, p.numEqu,2,p.numComp)];
519 const register double B_1_5 = B_p[INDEX4(k,1,m,5, p.numEqu,2,p.numComp)];
520 const register double B_0_6 = B_p[INDEX4(k,0,m,6, p.numEqu,2,p.numComp)];
521 const register double B_1_6 = B_p[INDEX4(k,1,m,6, p.numEqu,2,p.numComp)];
522 const register double B_0_7 = B_p[INDEX4(k,0,m,7, p.numEqu,2,p.numComp)];
523 const register double B_1_7 = B_p[INDEX4(k,1,m,7, p.numEqu,2,p.numComp)];
524 const register double B_0_8 = B_p[INDEX4(k,0,m,8, p.numEqu,2,p.numComp)];
525 const register double B_1_8 = B_p[INDEX4(k,1,m,8, p.numEqu,2,p.numComp)];
526 const register double tmp1_0 = B_1_6 + B_1_8;
527 const register double tmp2_0 = B_1_0 + B_1_2;
528 const register double tmp3_0 = B_0_2 + B_0_8;
529 const register double tmp5_0 = B_0_0 + B_0_6;
530 const register double tmp4_0 = B_0_1 + B_0_7;
531 const register double tmp0_0 = B_1_3 + B_1_5;
532 const register double tmp37_1 = B_1_0*w85;
533 const register double tmp89_1 = B_1_8*w85;
534 const register double tmp97_1 = B_0_6*w54;
535 const register double tmp83_1 = B_1_5*w90;
536 const register double tmp71_1 = tmp4_0*w93;
537 const register double tmp3_1 = B_0_1*w56;
538 const register double tmp4_1 = B_0_6*w64;
539 const register double tmp95_1 = B_0_8*w58;
540 const register double tmp38_1 = B_0_5*w59;
541 const register double tmp17_1 = B_0_5*w74;
542 const register double tmp64_1 = B_1_1*w84;
543 const register double tmp119_1 = B_0_6*w78;
544 const register double tmp102_1 = B_0_0*w68;
545 const register double tmp24_1 = B_1_1*w67;
546 const register double tmp74_1 = B_1_6*w85;
547 const register double tmp36_1 = B_0_2*w78;
548 const register double tmp77_1 = B_1_0*w55;
549 const register double tmp46_1 = B_1_8*w88;
550 const register double tmp2_1 = B_0_3*w59;
551 const register double tmp57_1 = B_1_3*w90;
552 const register double tmp52_1 = B_0_8*w81;
553 const register double tmp111_1 = B_0_7*w76;
554 const register double tmp92_1 = B_0_7*w56;
555 const register double tmp9_1 = tmp2_0*w55;
556 const register double tmp99_1 = B_1_8*w75;
557 const register double tmp0_1 = B_1_1*w57;
558 const register double tmp19_1 = tmp1_0*w55;
559 const register double tmp26_1 = B_1_4*w80;
560 const register double tmp33_1 = B_0_7*w83;
561 const register double tmp45_1 = B_0_2*w54;
562 const register double tmp44_1 = B_1_5*w87;
563 const register double tmp85_1 = B_1_8*w89;
564 const register double tmp79_1 = tmp5_0*w54;
565 const register double tmp15_1 = tmp3_0*w71;
566 const register double tmp112_1 = B_0_8*w71;
567 const register double tmp23_1 = B_0_3*w72;
568 const register double tmp109_1 = B_0_1*w83;
569 const register double tmp21_1 = B_0_4*w73;
570 const register double tmp54_1 = B_1_2*w82;
571 const register double tmp29_1 = tmp1_0*w82;
572 const register double tmp43_1 = B_1_6*w55;
573 const register double tmp47_1 = B_0_3*w63;
574 const register double tmp94_1 = B_0_2*w68;
575 const register double tmp116_1 = B_1_0*w65;
576 const register double tmp7_1 = B_0_4*w61;
577 const register double tmp56_1 = B_0_5*w72;
578 const register double tmp12_1 = B_0_0*w54;
579 const register double tmp30_1 = B_1_1*w77;
580 const register double tmp10_1 = B_1_4*w62;
581 const register double tmp13_1 = B_0_5*w63;
582 const register double tmp107_1 = B_0_0*w81;
583 const register double tmp66_1 = B_1_2*w92;
584 const register double tmp65_1 = tmp5_0*w71;
585 const register double tmp14_1 = B_1_7*w67;
586 const register double tmp106_1 = B_0_6*w71;
587 const register double tmp113_1 = B_1_2*w85;
588 const register double tmp63_1 = tmp2_0*w82;
589 const register double tmp93_1 = B_0_0*w64;
590 const register double tmp28_1 = B_1_7*w84;
591 const register double tmp75_1 = tmp3_0*w68;
592 const register double tmp80_1 = B_1_0*w92;
593 const register double tmp53_1 = B_1_8*w92;
594 const register double tmp49_1 = B_0_2*w71;
595 const register double tmp76_1 = B_1_8*w65;
596 const register double tmp40_1 = B_0_0*w58;
597 const register double tmp31_1 = B_0_1*w76;
598 const register double tmp86_1 = B_1_5*w86;
599 const register double tmp50_1 = B_0_6*w69;
600 const register double tmp22_1 = tmp5_0*w69;
601 const register double tmp90_1 = B_1_3*w87;
602 const register double tmp87_1 = B_1_6*w65;
603 const register double tmp88_1 = B_1_2*w55;
604 const register double tmp98_1 = B_1_6*w92;
605 const register double tmp6_1 = B_0_8*w68;
606 const register double tmp11_1 = B_0_7*w66;
607 const register double tmp58_1 = B_0_0*w78;
608 const register double tmp69_1 = tmp3_0*w69;
609 const register double tmp48_1 = B_0_6*w68;
610 const register double tmp68_1 = B_1_8*w82;
611 const register double tmp108_1 = B_0_2*w69;
612 const register double tmp96_1 = B_0_1*w66;
613 const register double tmp101_1 = B_1_0*w82;
614 const register double tmp34_1 = tmp0_0*w79;
615 const register double tmp70_1 = B_1_6*w89;
616 const register double tmp18_1 = tmp4_0*w70;
617 const register double tmp59_1 = B_1_6*w75;
618 const register double tmp72_1 = tmp5_0*w68;
619 const register double tmp51_1 = B_0_3*w74;
620 const register double tmp25_1 = B_0_0*w71;
621 const register double tmp118_1 = B_1_6*w88;
622 const register double tmp61_1 = tmp1_0*w75;
623 const register double tmp78_1 = B_1_2*w88;
624 const register double tmp16_1 = B_1_7*w57;
625 const register double tmp35_1 = tmp2_0*w75;
626 const register double tmp5_1 = B_0_2*w58;
627 const register double tmp105_1 = B_1_2*w89;
628 const register double tmp100_1 = B_0_2*w64;
629 const register double tmp84_1 = B_1_2*w75;
630 const register double tmp114_1 = B_0_2*w81;
631 const register double tmp62_1 = B_1_7*w77;
632 const register double tmp91_1 = B_1_0*w88;
633 const register double tmp110_1 = B_0_8*w78;
634 const register double tmp27_1 = B_0_6*w81;
635 const register double tmp32_1 = B_0_8*w69;
636 const register double tmp42_1 = B_1_2*w65;
637 const register double tmp8_1 = tmp1_0*w65;
638 const register double tmp1_1 = tmp0_0*w60;
639 const register double tmp104_1 = B_0_8*w54;
640 const register double tmp117_1 = B_0_0*w69;
641 const register double tmp67_1 = B_1_0*w75;
642 const register double tmp20_1 = tmp2_0*w65;
643 const register double tmp103_1 = B_0_6*w58;
644 const register double tmp115_1 = B_1_8*w55;
645 const register double tmp60_1 = B_1_0*w89;
646 const register double tmp41_1 = B_1_3*w86;
647 const register double tmp55_1 = B_1_5*w91;
648 const register double tmp82_1 = B_1_6*w82;
649 const register double tmp39_1 = B_0_8*w64;
650 const register double tmp73_1 = tmp3_0*w54;
651 const register double tmp81_1 = B_1_3*w91;
652 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
653 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp1_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1;
654 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp17_1 + tmp21_1 + tmp23_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp35_1 + tmp36_1;
655 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp14_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp3_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp7_1;
656 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp21_1 + tmp26_1 + tmp28_1 + tmp30_1 + tmp31_1 + tmp33_1 + tmp49_1 + tmp50_1 + tmp51_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1 + tmp60_1;
657 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp15_1 + tmp17_1 + tmp18_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp26_1 + tmp34_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1;
658 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp18_1 + tmp21_1 + tmp26_1 + tmp51_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp62_1 + tmp64_1 + tmp65_1 + tmp66_1 + tmp67_1 + tmp68_1 + tmp69_1 + tmp70_1;
659 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp13_1 + tmp26_1 + tmp2_1 + tmp34_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1 + tmp71_1 + tmp72_1 + tmp73_1 + tmp7_1;
660 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp16_1 + tmp24_1 + tmp38_1 + tmp41_1 + tmp44_1 + tmp47_1 + tmp71_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp78_1 + tmp79_1 + tmp7_1;
661 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp26_1 + tmp38_1 + tmp47_1 + tmp62_1 + tmp64_1 + tmp71_1 + tmp75_1 + tmp79_1 + tmp7_1 + tmp80_1 + tmp81_1 + tmp82_1 + tmp83_1 + tmp84_1 + tmp85_1;
662 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp16_1 + tmp18_1 + tmp21_1 + tmp24_1 + tmp51_1 + tmp56_1 + tmp65_1 + tmp69_1 + tmp86_1 + tmp87_1 + tmp88_1 + tmp89_1 + tmp90_1 + tmp91_1;
663 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp13_1 + tmp26_1 + tmp28_1 + tmp29_1 + tmp2_1 + tmp30_1 + tmp34_1 + tmp35_1 + tmp7_1 + tmp92_1 + tmp93_1 + tmp94_1 + tmp95_1 + tmp96_1 + tmp97_1;
664 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp100_1 + tmp101_1 + tmp102_1 + tmp103_1 + tmp104_1 + tmp105_1 + tmp26_1 + tmp28_1 + tmp30_1 + tmp38_1 + tmp47_1 + tmp7_1 + tmp81_1 + tmp83_1 + tmp92_1 + tmp96_1 + tmp98_1 + tmp99_1;
665 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp106_1 + tmp107_1 + tmp108_1 + tmp109_1 + tmp10_1 + tmp110_1 + tmp111_1 + tmp14_1 + tmp17_1 + tmp1_1 + tmp21_1 + tmp23_1 + tmp8_1 + tmp9_1;
666 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp13_1 + tmp16_1 + tmp19_1 + tmp1_1 + tmp20_1 + tmp24_1 + tmp2_1 + tmp71_1 + tmp72_1 + tmp73_1 + tmp7_1;
667 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp109_1 + tmp10_1 + tmp111_1 + tmp112_1 + tmp113_1 + tmp114_1 + tmp115_1 + tmp116_1 + tmp117_1 + tmp118_1 + tmp119_1 + tmp14_1 + tmp21_1 + tmp51_1 + tmp56_1 + tmp86_1 + tmp90_1;
668 }
669 }
670 } else { /* constant data */
671 for (k=0;k<p.numEqu;k++) {
672 for (m=0;m<p.numComp;m++) {
673 const register double B_0 = B_p[INDEX3(k,0,m, p.numEqu,2)];
674 const register double B_1 = B_p[INDEX3(k,1,m, p.numEqu,2)];
675 const register double tmp1_1 = B_1*w95;
676 const register double tmp2_1 = B_0*w96;
677 const register double tmp0_1 = B_0*w94;
678 const register double tmp6_1 = B_1*w100;
679 const register double tmp7_1 = B_0*w101;
680 const register double tmp4_1 = B_0*w97;
681 const register double tmp3_1 = B_1*w98;
682 const register double tmp5_1 = B_1*w99;
683 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
684 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp2_1;
685 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp3_1 + tmp4_1;
686 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp5_1;
687 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1;
688 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
689 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp6_1;
690 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp3_1 + tmp7_1;
691 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp5_1 + tmp7_1;
692 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp6_1 + tmp7_1;
693 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp5_1;
694 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp3_1;
695 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp6_1;
696 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp4_1;
697 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp7_1;
698 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1;
699 }
700 }
701 }
702 }
703 /**************************************************************/
704 /* process C: */
705 /**************************************************************/
706 if (NULL!=C_p) {
707 add_EM_S=TRUE;
708 if (extendedC) {
709 for (k=0;k<p.numEqu;k++) {
710 for (m=0;m<p.numComp;m++) {
711 const register double C_0_0 = C_p[INDEX4(k,m,0, 0, p.numEqu,p.numComp,2)];
712 const register double C_1_0 = C_p[INDEX4(k,m,1, 0, p.numEqu,p.numComp,2)];
713 const register double C_0_1 = C_p[INDEX4(k,m,0, 1, p.numEqu,p.numComp,2)];
714 const register double C_1_1 = C_p[INDEX4(k,m,1, 1, p.numEqu,p.numComp,2)];
715 const register double C_0_2 = C_p[INDEX4(k,m,0, 2, p.numEqu,p.numComp,2)];
716 const register double C_1_2 = C_p[INDEX4(k,m,1, 2, p.numEqu,p.numComp,2)];
717 const register double C_0_3 = C_p[INDEX4(k,m,0, 3, p.numEqu,p.numComp,2)];
718 const register double C_1_3 = C_p[INDEX4(k,m,1, 3, p.numEqu,p.numComp,2)];
719 const register double C_0_4 = C_p[INDEX4(k,m,0, 4, p.numEqu,p.numComp,2)];
720 const register double C_1_4 = C_p[INDEX4(k,m,1, 4, p.numEqu,p.numComp,2)];
721 const register double C_0_5 = C_p[INDEX4(k,m,0, 5, p.numEqu,p.numComp,2)];
722 const register double C_1_5 = C_p[INDEX4(k,m,1, 5, p.numEqu,p.numComp,2)];
723 const register double C_0_6 = C_p[INDEX4(k,m,0, 6, p.numEqu,p.numComp,2)];
724 const register double C_1_6 = C_p[INDEX4(k,m,1, 6, p.numEqu,p.numComp,2)];
725 const register double C_0_7 = C_p[INDEX4(k,m,0, 7, p.numEqu,p.numComp,2)];
726 const register double C_1_7 = C_p[INDEX4(k,m,1, 7, p.numEqu,p.numComp,2)];
727 const register double C_0_8 = C_p[INDEX4(k,m,0, 8, p.numEqu,p.numComp,2)];
728 const register double C_1_8 = C_p[INDEX4(k,m,1, 8, p.numEqu,p.numComp,2)];
729 const register double tmp5_0 = C_0_2 + C_0_8;
730 const register double tmp2_0 = C_1_0 + C_1_2;
731 const register double tmp0_0 = C_1_3 + C_1_5;
732 const register double tmp3_0 = C_0_1 + C_0_7;
733 const register double tmp4_0 = C_0_0 + C_0_6;
734 const register double tmp1_0 = C_1_6 + C_1_8;
735 const register double tmp68_1 = C_1_1*w67;
736 const register double tmp110_1 = C_0_8*w68;
737 const register double tmp85_1 = tmp4_0*w54;
738 const register double tmp25_1 = C_0_5*w63;
739 const register double tmp43_1 = C_1_2*w65;
740 const register double tmp19_1 = C_0_4*w61;
741 const register double tmp41_1 = C_0_0*w58;
742 const register double tmp9_1 = C_0_2*w69;
743 const register double tmp103_1 = C_1_0*w82;
744 const register double tmp47_1 = C_0_2*w54;
745 const register double tmp51_1 = C_0_2*w71;
746 const register double tmp95_1 = C_1_6*w89;
747 const register double tmp72_1 = C_1_6*w65;
748 const register double tmp7_1 = tmp2_0*w55;
749 const register double tmp4_1 = C_0_0*w81;
750 const register double tmp30_1 = C_0_2*w68;
751 const register double tmp79_1 = tmp4_0*w69;
752 const register double tmp82_1 = C_1_3*w91;
753 const register double tmp63_1 = C_1_6*w75;
754 const register double tmp71_1 = C_1_5*w86;
755 const register double tmp67_1 = tmp1_0*w55;
756 const register double tmp20_1 = tmp4_0*w68;
757 const register double tmp24_1 = tmp5_0*w54;
758 const register double tmp118_1 = C_1_6*w88;
759 const register double tmp97_1 = C_0_6*w81;
760 const register double tmp96_1 = C_0_0*w71;
761 const register double tmp70_1 = tmp3_0*w70;
762 const register double tmp34_1 = C_0_1*w66;
763 const register double tmp113_1 = C_1_2*w85;
764 const register double tmp99_1 = C_0_2*w78;
765 const register double tmp54_1 = C_0_8*w81;
766 const register double tmp42_1 = C_1_3*w86;
767 const register double tmp87_1 = C_1_8*w89;
768 const register double tmp61_1 = C_1_3*w90;
769 const register double tmp48_1 = C_1_8*w88;
770 const register double tmp55_1 = C_1_8*w92;
771 const register double tmp0_1 = C_1_1*w57;
772 const register double tmp65_1 = C_1_7*w57;
773 const register double tmp94_1 = C_1_8*w82;
774 const register double tmp116_1 = C_1_0*w65;
775 const register double tmp1_1 = C_0_6*w71;
776 const register double tmp101_1 = C_1_8*w75;
777 const register double tmp32_1 = tmp1_0*w82;
778 const register double tmp2_1 = C_0_5*w74;
779 const register double tmp108_1 = C_0_6*w64;
780 const register double tmp8_1 = C_1_4*w62;
781 const register double tmp90_1 = C_1_0*w55;
782 const register double tmp14_1 = C_0_7*w76;
783 const register double tmp22_1 = C_1_7*w77;
784 const register double tmp40_1 = C_0_8*w64;
785 const register double tmp5_1 = C_0_4*w73;
786 const register double tmp74_1 = C_1_8*w85;
787 const register double tmp78_1 = tmp5_0*w71;
788 const register double tmp76_1 = C_1_3*w87;
789 const register double tmp57_1 = C_0_7*w83;
790 const register double tmp21_1 = tmp2_0*w82;
791 const register double tmp52_1 = C_0_6*w69;
792 const register double tmp39_1 = C_0_1*w56;
793 const register double tmp88_1 = C_1_6*w85;
794 const register double tmp31_1 = C_0_8*w58;
795 const register double tmp38_1 = C_0_5*w59;
796 const register double tmp10_1 = C_0_3*w72;
797 const register double tmp102_1 = C_0_2*w64;
798 const register double tmp80_1 = C_1_0*w92;
799 const register double tmp49_1 = C_0_3*w63;
800 const register double tmp75_1 = tmp5_0*w69;
801 const register double tmp77_1 = C_1_0*w88;
802 const register double tmp104_1 = C_0_0*w68;
803 const register double tmp60_1 = C_0_5*w72;
804 const register double tmp56_1 = C_1_2*w82;
805 const register double tmp69_1 = tmp4_0*w71;
806 const register double tmp91_1 = C_1_2*w88;
807 const register double tmp93_1 = C_1_0*w75;
808 const register double tmp107_1 = C_1_2*w89;
809 const register double tmp106_1 = C_0_8*w54;
810 const register double tmp13_1 = C_0_8*w78;
811 const register double tmp64_1 = C_1_0*w89;
812 const register double tmp33_1 = C_1_1*w77;
813 const register double tmp115_1 = C_1_8*w55;
814 const register double tmp11_1 = C_0_1*w83;
815 const register double tmp86_1 = C_1_2*w75;
816 const register double tmp83_1 = C_1_6*w82;
817 const register double tmp109_1 = C_0_2*w58;
818 const register double tmp98_1 = C_0_8*w69;
819 const register double tmp23_1 = tmp0_0*w79;
820 const register double tmp114_1 = C_0_2*w81;
821 const register double tmp81_1 = tmp5_0*w68;
822 const register double tmp53_1 = C_0_3*w74;
823 const register double tmp92_1 = C_1_2*w92;
824 const register double tmp119_1 = C_0_6*w78;
825 const register double tmp50_1 = C_0_6*w68;
826 const register double tmp18_1 = tmp1_0*w75;
827 const register double tmp15_1 = tmp3_0*w93;
828 const register double tmp73_1 = C_1_2*w55;
829 const register double tmp36_1 = tmp2_0*w75;
830 const register double tmp17_1 = C_0_3*w59;
831 const register double tmp45_1 = C_1_5*w87;
832 const register double tmp27_1 = C_0_7*w56;
833 const register double tmp84_1 = C_1_5*w90;
834 const register double tmp59_1 = C_1_5*w91;
835 const register double tmp105_1 = C_0_6*w58;
836 const register double tmp44_1 = C_1_6*w55;
837 const register double tmp89_1 = C_1_8*w65;
838 const register double tmp62_1 = C_0_0*w78;
839 const register double tmp12_1 = C_1_7*w67;
840 const register double tmp100_1 = C_1_6*w92;
841 const register double tmp46_1 = C_0_7*w66;
842 const register double tmp117_1 = C_0_0*w69;
843 const register double tmp6_1 = tmp1_0*w65;
844 const register double tmp3_1 = tmp0_0*w60;
845 const register double tmp37_1 = C_1_0*w85;
846 const register double tmp111_1 = C_0_0*w54;
847 const register double tmp66_1 = tmp2_0*w65;
848 const register double tmp35_1 = C_0_6*w54;
849 const register double tmp112_1 = C_0_8*w71;
850 const register double tmp29_1 = C_0_0*w64;
851 const register double tmp28_1 = C_1_7*w84;
852 const register double tmp16_1 = C_1_4*w80;
853 const register double tmp26_1 = C_1_1*w84;
854 const register double tmp58_1 = C_0_1*w76;
855 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
856 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp26_1;
857 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp16_1 + tmp17_1 + tmp19_1 + tmp23_1 + tmp25_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp35_1 + tmp36_1;
858 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp12_1 + tmp19_1 + tmp37_1 + tmp38_1 + tmp39_1 + tmp40_1 + tmp41_1 + tmp42_1 + tmp43_1 + tmp44_1 + tmp45_1 + tmp46_1 + tmp47_1 + tmp48_1 + tmp49_1 + tmp50_1 + tmp8_1;
859 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp16_1 + tmp28_1 + tmp33_1 + tmp51_1 + tmp52_1 + tmp53_1 + tmp54_1 + tmp55_1 + tmp56_1 + tmp57_1 + tmp58_1 + tmp59_1 + tmp5_1 + tmp60_1 + tmp61_1 + tmp62_1 + tmp63_1 + tmp64_1;
860 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp15_1 + tmp17_1 + tmp19_1 + tmp20_1 + tmp24_1 + tmp25_1 + tmp3_1 + tmp65_1 + tmp66_1 + tmp67_1 + tmp68_1 + tmp8_1;
861 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp53_1 + tmp5_1 + tmp60_1 + tmp65_1 + tmp68_1 + tmp69_1 + tmp70_1 + tmp71_1 + tmp72_1 + tmp73_1 + tmp74_1 + tmp75_1 + tmp76_1 + tmp77_1 + tmp8_1;
862 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp2_1 + tmp3_1 + tmp5_1 + tmp65_1 + tmp66_1 + tmp67_1 + tmp68_1 + tmp70_1 + tmp78_1 + tmp79_1 + tmp8_1;
863 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp15_1 + tmp16_1 + tmp19_1 + tmp22_1 + tmp26_1 + tmp38_1 + tmp49_1 + tmp80_1 + tmp81_1 + tmp82_1 + tmp83_1 + tmp84_1 + tmp85_1 + tmp86_1 + tmp87_1;
864 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp15_1 + tmp19_1 + tmp38_1 + tmp42_1 + tmp45_1 + tmp49_1 + tmp65_1 + tmp68_1 + tmp81_1 + tmp85_1 + tmp88_1 + tmp89_1 + tmp8_1 + tmp90_1 + tmp91_1;
865 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp16_1 + tmp22_1 + tmp26_1 + tmp53_1 + tmp59_1 + tmp5_1 + tmp60_1 + tmp61_1 + tmp69_1 + tmp70_1 + tmp75_1 + tmp92_1 + tmp93_1 + tmp94_1 + tmp95_1;
866 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp16_1 + tmp23_1 + tmp28_1 + tmp2_1 + tmp32_1 + tmp33_1 + tmp36_1 + tmp57_1 + tmp58_1 + tmp5_1 + tmp96_1 + tmp97_1 + tmp98_1 + tmp99_1;
867 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp100_1 + tmp101_1 + tmp102_1 + tmp103_1 + tmp104_1 + tmp105_1 + tmp106_1 + tmp107_1 + tmp16_1 + tmp19_1 + tmp27_1 + tmp28_1 + tmp33_1 + tmp34_1 + tmp38_1 + tmp49_1 + tmp82_1 + tmp84_1;
868 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp108_1 + tmp109_1 + tmp110_1 + tmp111_1 + tmp12_1 + tmp17_1 + tmp19_1 + tmp25_1 + tmp39_1 + tmp3_1 + tmp46_1 + tmp6_1 + tmp7_1 + tmp8_1;
869 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp16_1 + tmp18_1 + tmp21_1 + tmp22_1 + tmp23_1 + tmp26_1 + tmp2_1 + tmp5_1 + tmp70_1 + tmp78_1 + tmp79_1;
870 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp112_1 + tmp113_1 + tmp114_1 + tmp115_1 + tmp116_1 + tmp117_1 + tmp118_1 + tmp119_1 + tmp11_1 + tmp12_1 + tmp14_1 + tmp53_1 + tmp5_1 + tmp60_1 + tmp71_1 + tmp76_1 + tmp8_1;
871 }
872 }
873 } else { /* constant data */
874 for (k=0;k<p.numEqu;k++) {
875 for (m=0;m<p.numComp;m++) {
876 const register double C_0 = C_p[INDEX3(k,m,0, p.numEqu,p.numComp)];
877 const register double C_1 = C_p[INDEX3(k,m,1, p.numEqu,p.numComp)];
878 const register double tmp1_1 = C_1*w95;
879 const register double tmp3_1 = C_0*w101;
880 const register double tmp0_1 = C_0*w97;
881 const register double tmp6_1 = C_1*w100;
882 const register double tmp7_1 = C_0*w96;
883 const register double tmp4_1 = C_0*w94;
884 const register double tmp2_1 = C_1*w98;
885 const register double tmp5_1 = C_1*w99;
886 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1;
887 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp3_1;
888 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp4_1;
889 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp5_1;
890 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp6_1;
891 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp3_1;
892 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp5_1 + tmp7_1;
893 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp7_1;
894 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp3_1 + tmp6_1;
895 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp3_1 + tmp5_1;
896 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp6_1 + tmp7_1;
897 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp2_1;
898 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1;
899 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp1_1 + tmp4_1;
900 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp2_1 + tmp7_1;
901 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp5_1;
902 }
903 }
904 }
905 }
906 /**************************************************************/
907 /* process D: */
908 /**************************************************************/
909 if (NULL!=D_p) {
910 add_EM_S=TRUE;
911 if (extendedD) {
912 for (k=0;k<p.numEqu;k++) {
913 for (m=0;m<p.numComp;m++) {
914 const register double D_0 = D_p[INDEX3(k,m,0, p.numEqu,p.numComp)];
915 const register double D_1 = D_p[INDEX3(k,m,1, p.numEqu,p.numComp)];
916 const register double D_2 = D_p[INDEX3(k,m,2, p.numEqu,p.numComp)];
917 const register double D_3 = D_p[INDEX3(k,m,3, p.numEqu,p.numComp)];
918 const register double D_4 = D_p[INDEX3(k,m,4, p.numEqu,p.numComp)];
919 const register double D_5 = D_p[INDEX3(k,m,5, p.numEqu,p.numComp)];
920 const register double D_6 = D_p[INDEX3(k,m,6, p.numEqu,p.numComp)];
921 const register double D_7 = D_p[INDEX3(k,m,7, p.numEqu,p.numComp)];
922 const register double D_8 = D_p[INDEX3(k,m,8, p.numEqu,p.numComp)];
923 const register double tmp12_0 = D_1 + D_5;
924 const register double tmp7_0 = D_1 + D_3;
925 const register double tmp3_0 = D_1 + D_3 + D_5 + D_7;
926 const register double tmp11_0 = D_0 + D_8;
927 const register double tmp2_0 = D_0 + D_2;
928 const register double tmp9_0 = D_1 + D_7;
929 const register double tmp0_0 = D_6 + D_8;
930 const register double tmp8_0 = D_0 + D_6;
931 const register double tmp10_0 = D_2 + D_8;
932 const register double tmp1_0 = D_3 + D_5;
933 const register double tmp6_0 = D_5 + D_7;
934 const register double tmp5_0 = D_2 + D_6;
935 const register double tmp4_0 = D_0 + D_2 + D_6 + D_8;
936 const register double tmp13_0 = D_3 + D_7;
937 const register double tmp2_1 = D_7*w107;
938 const register double tmp18_1 = tmp7_0*w107;
939 const register double tmp5_1 = tmp2_0*w102;
940 const register double tmp25_1 = tmp10_0*w102;
941 const register double tmp32_1 = tmp12_0*w107;
942 const register double tmp30_1 = tmp11_0*w108;
943 const register double tmp35_1 = D_6*w110;
944 const register double tmp31_1 = D_2*w110;
945 const register double tmp22_1 = D_5*w103;
946 const register double tmp24_1 = tmp9_0*w104;
947 const register double tmp29_1 = tmp8_0*w102;
948 const register double tmp38_1 = tmp12_0*w103;
949 const register double tmp9_1 = D_7*w103;
950 const register double tmp23_1 = D_3*w107;
951 const register double tmp1_1 = D_1*w103;
952 const register double tmp27_1 = D_3*w103;
953 const register double tmp34_1 = tmp13_0*w103;
954 const register double tmp15_1 = D_0*w109;
955 const register double tmp36_1 = tmp13_0*w107;
956 const register double tmp16_1 = tmp7_0*w103;
957 const register double tmp19_1 = D_8*w109;
958 const register double tmp26_1 = tmp10_0*w106;
959 const register double tmp33_1 = D_6*w109;
960 const register double tmp12_1 = tmp5_0*w108;
961 const register double tmp8_1 = tmp2_0*w106;
962 const register double tmp13_1 = D_8*w110;
963 const register double tmp6_1 = tmp3_0*w104;
964 const register double tmp20_1 = D_0*w110;
965 const register double tmp28_1 = D_5*w107;
966 const register double tmp0_1 = tmp0_0*w106;
967 const register double tmp37_1 = D_2*w109;
968 const register double tmp17_1 = tmp6_0*w103;
969 const register double tmp14_1 = tmp6_0*w107;
970 const register double tmp21_1 = tmp8_0*w106;
971 const register double tmp3_1 = tmp1_0*w104;
972 const register double tmp7_1 = tmp4_0*w108;
973 const register double tmp11_1 = tmp0_0*w102;
974 const register double tmp4_1 = D_4*w105;
975 const register double tmp10_1 = D_1*w107;
976 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
977 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1 + tmp7_1;
978 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp11_1 + tmp3_1 + tmp4_1 + tmp8_1 + tmp9_1;
979 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp4_1;
980 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp12_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp4_1;
981 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1 + tmp7_1;
982 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp4_1;
983 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1 + tmp7_1;
984 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp24_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp4_1;
985 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp24_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp4_1;
986 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp21_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp4_1;
987 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp10_1 + tmp11_1 + tmp3_1 + tmp4_1 + tmp8_1 + tmp9_1;
988 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp30_1 + tmp31_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp4_1;
989 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
990 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp4_1 + tmp6_1 + tmp7_1;
991 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp30_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp4_1;
992 }
993 }
994 } else { /* constant data */
995 for (k=0;k<p.numEqu;k++) {
996 for (m=0;m<p.numComp;m++) {
997 const register double D_0 = D_p[INDEX2(k,m, p.numEqu)];
998 const register double tmp0_1 = D_0*w111;
999 const register double tmp2_1 = D_0*w113;
1000 const register double tmp1_1 = D_0*w112;
1001 EM_S[INDEX4(k,m,0,1,p.numEqu,p.numComp,4)]+=tmp0_1;
1002 EM_S[INDEX4(k,m,1,2,p.numEqu,p.numComp,4)]+=tmp1_1;
1003 EM_S[INDEX4(k,m,3,2,p.numEqu,p.numComp,4)]+=tmp0_1;
1004 EM_S[INDEX4(k,m,0,0,p.numEqu,p.numComp,4)]+=tmp2_1;
1005 EM_S[INDEX4(k,m,3,3,p.numEqu,p.numComp,4)]+=tmp2_1;
1006 EM_S[INDEX4(k,m,3,0,p.numEqu,p.numComp,4)]+=tmp1_1;
1007 EM_S[INDEX4(k,m,3,1,p.numEqu,p.numComp,4)]+=tmp0_1;
1008 EM_S[INDEX4(k,m,2,1,p.numEqu,p.numComp,4)]+=tmp1_1;
1009 EM_S[INDEX4(k,m,0,2,p.numEqu,p.numComp,4)]+=tmp0_1;
1010 EM_S[INDEX4(k,m,2,0,p.numEqu,p.numComp,4)]+=tmp0_1;
1011 EM_S[INDEX4(k,m,1,3,p.numEqu,p.numComp,4)]+=tmp0_1;
1012 EM_S[INDEX4(k,m,2,3,p.numEqu,p.numComp,4)]+=tmp0_1;
1013 EM_S[INDEX4(k,m,2,2,p.numEqu,p.numComp,4)]+=tmp2_1;
1014 EM_S[INDEX4(k,m,1,0,p.numEqu,p.numComp,4)]+=tmp0_1;
1015 EM_S[INDEX4(k,m,0,3,p.numEqu,p.numComp,4)]+=tmp1_1;
1016 EM_S[INDEX4(k,m,1,1,p.numEqu,p.numComp,4)]+=tmp2_1;
1017 }
1018 }
1019 }
1020 }
1021 /**************************************************************/
1022 /* process X: */
1023 /**************************************************************/
1024 if (NULL!=X_p) {
1025 add_EM_F=TRUE;
1026 if (extendedX) {
1027 for (k=0;k<p.numEqu;k++) {
1028 const register double X_0_0 = X_p[INDEX3(k,0, 0, p.numEqu,2)];
1029 const register double X_1_0 = X_p[INDEX3(k,1, 0, p.numEqu,2)];
1030 const register double X_0_1 = X_p[INDEX3(k,0, 1, p.numEqu,2)];
1031 const register double X_1_1 = X_p[INDEX3(k,1, 1, p.numEqu,2)];
1032 const register double X_0_2 = X_p[INDEX3(k,0, 2, p.numEqu,2)];
1033 const register double X_1_2 = X_p[INDEX3(k,1, 2, p.numEqu,2)];
1034 const register double X_0_3 = X_p[INDEX3(k,0, 3, p.numEqu,2)];
1035 const register double X_1_3 = X_p[INDEX3(k,1, 3, p.numEqu,2)];
1036 const register double X_0_4 = X_p[INDEX3(k,0, 4, p.numEqu,2)];
1037 const register double X_1_4 = X_p[INDEX3(k,1, 4, p.numEqu,2)];
1038 const register double X_0_5 = X_p[INDEX3(k,0, 5, p.numEqu,2)];
1039 const register double X_1_5 = X_p[INDEX3(k,1, 5, p.numEqu,2)];
1040 const register double X_0_6 = X_p[INDEX3(k,0, 6, p.numEqu,2)];
1041 const register double X_1_6 = X_p[INDEX3(k,1, 6, p.numEqu,2)];
1042 const register double X_0_7 = X_p[INDEX3(k,0, 7, p.numEqu,2)];
1043 const register double X_1_7 = X_p[INDEX3(k,1, 7, p.numEqu,2)];
1044 const register double X_0_8 = X_p[INDEX3(k,0, 8, p.numEqu,2)];
1045 const register double X_1_8 = X_p[INDEX3(k,1, 8, p.numEqu,2)];
1046 const register double tmp0_0 = X_1_0 + X_1_6;
1047 const register double tmp4_0 = X_0_3 + X_0_5;
1048 const register double tmp2_0 = X_0_0 + X_0_2;
1049 const register double tmp5_0 = X_1_2 + X_1_8;
1050 const register double tmp1_0 = X_0_6 + X_0_8;
1051 const register double tmp3_0 = X_1_1 + X_1_7;
1052 const register double tmp24_1 = X_1_3*w135;
1053 const register double tmp14_1 = tmp2_0*w126;
1054 const register double tmp8_1 = X_0_1*w116;
1055 const register double tmp22_1 = tmp2_0*w124;
1056 const register double tmp34_1 = X_1_5*w135;
1057 const register double tmp15_1 = X_0_4*w129;
1058 const register double tmp20_1 = X_1_5*w120;
1059 const register double tmp36_1 = tmp0_0*w134;
1060 const register double tmp39_1 = tmp5_0*w132;
1061 const register double tmp0_1 = X_0_4*w121;
1062 const register double tmp5_1 = X_1_5*w123;
1063 const register double tmp38_1 = X_1_3*w137;
1064 const register double tmp3_1 = X_1_4*w122;
1065 const register double tmp2_1 = tmp1_0*w124;
1066 const register double tmp11_1 = X_0_7*w125;
1067 const register double tmp19_1 = X_0_7*w131;
1068 const register double tmp1_1 = tmp0_0*w115;
1069 const register double tmp10_1 = tmp5_0*w118;
1070 const register double tmp12_1 = tmp5_0*w115;
1071 const register double tmp16_1 = X_0_1*w127;
1072 const register double tmp23_1 = X_1_4*w136;
1073 const register double tmp31_1 = X_0_1*w125;
1074 const register double tmp17_1 = tmp4_0*w128;
1075 const register double tmp37_1 = tmp1_0*w126;
1076 const register double tmp30_1 = tmp0_0*w132;
1077 const register double tmp26_1 = tmp5_0*w134;
1078 const register double tmp29_1 = X_1_5*w137;
1079 const register double tmp7_1 = tmp4_0*w119;
1080 const register double tmp4_1 = tmp2_0*w114;
1081 const register double tmp18_1 = X_1_3*w123;
1082 const register double tmp33_1 = X_0_1*w131;
1083 const register double tmp9_1 = X_1_3*w120;
1084 const register double tmp32_1 = X_0_7*w127;
1085 const register double tmp13_1 = tmp1_0*w130;
1086 const register double tmp35_1 = tmp2_0*w130;
1087 const register double tmp6_1 = tmp3_0*w117;
1088 const register double tmp28_1 = X_0_7*w116;
1089 const register double tmp21_1 = tmp0_0*w118;
1090 const register double tmp27_1 = tmp3_0*w133;
1091 const register double tmp25_1 = tmp1_0*w114;
1092 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1 + tmp10_1 + tmp11_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1093 EM_F[INDEX2(k,1,p.numEqu)]+=tmp12_1 + tmp13_1 + tmp14_1 + tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp19_1 + tmp20_1 + tmp21_1 + tmp3_1 + tmp6_1;
1094 EM_F[INDEX2(k,2,p.numEqu)]+=tmp0_1 + tmp22_1 + tmp23_1 + tmp24_1 + tmp25_1 + tmp26_1 + tmp27_1 + tmp28_1 + tmp29_1 + tmp30_1 + tmp31_1 + tmp7_1;
1095 EM_F[INDEX2(k,3,p.numEqu)]+=tmp15_1 + tmp17_1 + tmp23_1 + tmp27_1 + tmp32_1 + tmp33_1 + tmp34_1 + tmp35_1 + tmp36_1 + tmp37_1 + tmp38_1 + tmp39_1;
1096 }
1097 } else { /* constant data */
1098 for (k=0;k<p.numEqu;k++) {
1099 const register double X_0 = X_p[INDEX2(k,0, p.numEqu)];
1100 const register double X_1 = X_p[INDEX2(k,1, p.numEqu)];
1101 const register double tmp2_1 = X_0*w140;
1102 const register double tmp3_1 = X_1*w141;
1103 const register double tmp1_1 = X_1*w139;
1104 const register double tmp0_1 = X_0*w138;
1105 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1 + tmp1_1;
1106 EM_F[INDEX2(k,1,p.numEqu)]+=tmp1_1 + tmp2_1;
1107 EM_F[INDEX2(k,2,p.numEqu)]+=tmp0_1 + tmp3_1;
1108 EM_F[INDEX2(k,3,p.numEqu)]+=tmp2_1 + tmp3_1;
1109 }
1110 }
1111 }
1112 /**************************************************************/
1113 /* process Y: */
1114 /**************************************************************/
1115 if (NULL!=Y_p) {
1116 add_EM_F=TRUE;
1117 if (extendedY) {
1118 for (k=0;k<p.numEqu;k++) {
1119 const register double Y_0 = Y_p[INDEX3(k,0, p.numEqu)];
1120 const register double Y_1 = Y_p[INDEX3(k,1, p.numEqu)];
1121 const register double Y_2 = Y_p[INDEX3(k,2, p.numEqu)];
1122 const register double Y_3 = Y_p[INDEX3(k,3, p.numEqu)];
1123 const register double Y_4 = Y_p[INDEX3(k,4, p.numEqu)];
1124 const register double Y_5 = Y_p[INDEX3(k,5, p.numEqu)];
1125 const register double Y_6 = Y_p[INDEX3(k,6, p.numEqu)];
1126 const register double Y_7 = Y_p[INDEX3(k,7, p.numEqu)];
1127 const register double Y_8 = Y_p[INDEX3(k,8, p.numEqu)];
1128 const register double tmp5_0 = Y_0 + Y_8;
1129 const register double tmp0_0 = Y_5 + Y_7;
1130 const register double tmp1_0 = Y_1 + Y_3;
1131 const register double tmp4_0 = Y_1 + Y_5;
1132 const register double tmp3_0 = Y_3 + Y_7;
1133 const register double tmp2_0 = Y_2 + Y_6;
1134 const register double tmp14_1 = tmp3_0*w143;
1135 const register double tmp11_1 = tmp4_0*w146;
1136 const register double tmp6_1 = tmp3_0*w146;
1137 const register double tmp17_1 = Y_0*w147;
1138 const register double tmp16_1 = Y_8*w142;
1139 const register double tmp13_1 = Y_2*w147;
1140 const register double tmp18_1 = tmp0_0*w143;
1141 const register double tmp10_1 = tmp5_0*w144;
1142 const register double tmp3_1 = tmp1_0*w143;
1143 const register double tmp12_1 = Y_6*w142;
1144 const register double tmp0_1 = tmp0_0*w146;
1145 const register double tmp9_1 = tmp4_0*w143;
1146 const register double tmp4_1 = tmp2_0*w144;
1147 const register double tmp2_1 = Y_8*w147;
1148 const register double tmp15_1 = tmp1_0*w146;
1149 const register double tmp8_1 = Y_6*w147;
1150 const register double tmp7_1 = Y_2*w142;
1151 const register double tmp5_1 = Y_4*w145;
1152 const register double tmp1_1 = Y_0*w142;
1153 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1 + tmp1_1 + tmp2_1 + tmp3_1 + tmp4_1 + tmp5_1;
1154 EM_F[INDEX2(k,1,p.numEqu)]+=tmp10_1 + tmp5_1 + tmp6_1 + tmp7_1 + tmp8_1 + tmp9_1;
1155 EM_F[INDEX2(k,2,p.numEqu)]+=tmp10_1 + tmp11_1 + tmp12_1 + tmp13_1 + tmp14_1 + tmp5_1;
1156 EM_F[INDEX2(k,3,p.numEqu)]+=tmp15_1 + tmp16_1 + tmp17_1 + tmp18_1 + tmp4_1 + tmp5_1;
1157 }
1158 } else { /* constant data */
1159 for (k=0;k<p.numEqu;k++) {
1160 const register double Y_0 = Y_p[k];
1161 const register double tmp0_1 = Y_0*w148;
1162 EM_F[INDEX2(k,0,p.numEqu)]+=tmp0_1;
1163 EM_F[INDEX2(k,1,p.numEqu)]+=tmp0_1;
1164 EM_F[INDEX2(k,2,p.numEqu)]+=tmp0_1;
1165 EM_F[INDEX2(k,3,p.numEqu)]+=tmp0_1;
1166 }
1167 }
1168 }
1169 /* GENERATOR SNIP BOTTOM */
1170 /* add to matrix (if add_EM_S) and RHS (if add_EM_F)*/
1171 } /* end k0 */
1172 } /* end k1 */
1173 } /* end k2 */
1174 } /* end coloring */
1175 } /* end of pointer check */
1176 THREAD_MEMFREE(EM_S);
1177 THREAD_MEMFREE(EM_F);
1178 } /* end parallel region */
1179 } /* end of presettings */
1180 }

  ViewVC Help
Powered by ViewVC 1.1.26