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

Annotation of /trunk/ripley/src/Assemble_PDE_System_3D.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3641 - (hide annotations)
Thu Oct 27 02:16:12 2011 UTC (7 years, 10 months ago) by gross
File MIME type: text/plain
File size: 98657 byte(s)
more work on the assemblage
1 gross 3641 /*******************************************************
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 has p.numComp components in a 3D 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_3D(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 3
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