31
32
33
34#include "implicit_f.inc"
35
36
37
38 INTEGER N0,N1,N2,NC, LLL(*),JLL(*),SLL(*),IADLL(*)
40 . x(3,*),xll(*),sk(9),l1(3),l2(3),
alpha
41
42
43
44 INTEGER IK,IC,IAD
46 . xt1(3),yt1(3),zt1(3),xt2(3),yt2(3),zt2(3),
47 . x0,y0,z0,x1,y1,z1,x2,y2,z2,xr,yr,zr,xm,ym,zm,
norm
48
49
50
51
52
53
54
55
56
57
58
59
60
61 xt1(1) = sk(1)*l1(1)+sk(4)*l1(2)+sk(7)*l1(3)
62 xt1(2) = sk(2)*l1(1)+sk(5)*l1(2)+sk(8)*l1(3)
63 xt1(3) = sk(3)*l1(1)+sk(6)*l1(2)+sk(9)*l1(3)
64 norm = one/sqrt(xt1(1)*xt1(1)+xt1(2)*xt1(2)+xt1(3)*xt1(3))
68 IF (abs(xt1(1))>zep99) THEN
69 zt1(1) =-xt1(3)
70 zt1(2) = zero
71 zt1(3) = xt1(1)
72 ELSE
73 zt1(1) = zero
74 zt1(2) =-xt1(3)
75 zt1(3) = xt1(2)
76 ENDIF
77 norm = one/sqrt(zt1(1)*zt1(1)+zt1(2)*zt1(2)+zt1(3)*zt1(3))
81 yt1(1) = zt1(2)*xt1(3)-zt1(3)*xt1(2)
82 yt1(2) = zt1(3)*xt1(1)-zt1(1)*xt1(3)
83 yt1(3) = zt1(1)*xt1(2)-zt1(2)*xt1(1)
84
85 xt2(1) = sk(1)*l2(1)+sk(4)*l2(2)+sk(7)*l2(3)
86 xt2(2) = sk(2)*l2(1)+sk(5)*l2(2)+sk(8)*l2(3)
87 xt2(3) = sk(3)*l2(1)+sk(6)*l2(2)+sk(9)*l2(3)
88 norm = one/sqrt(xt2(1)*xt2(1)+xt2(2)*xt2(2)+xt2(3)*xt2(3))
92 IF (abs(xt2(1))>zep99) THEN
93 zt2(1) =-xt2(3)
94 zt2(2) = zero
95 zt2(3) = xt2(1)
96 ELSE
97 zt2(1) = zero
98 zt2(2) =-xt2(3)
99 zt2(3) = xt2(2)
100 ENDIF
101 norm = one/sqrt(zt1(1)*zt1(1)+zt1(2)*zt1(2)+zt1(3)*zt1(3))
105 yt2(1) = zt2(2)*xt2(3)-zt2(3)*xt2(2)
106 yt2(2) = zt2(3)*xt2(1)-zt2(1)*xt2(3)
107 yt2(3) = zt2(1)*xt2(2)-zt2(2)*xt2(1)
108
109 x1=x(1,n1)
110 y1=x(2,n1)
111 z1=x(3,n1)
112 x2=x(1,n2)
113 y2=x(2,n2)
114 z2=x(3,n2)
115 x0=half*(x1+x2)
116 y0=half*(y1+y2)
117 z0=half*(z1+z2)
118 xr = x1 - x0
119 yr = y1 - y0
120 zr = z1 - z0
121 xm = x(1,n0) - x0
122 ym = x(2,n0) - y0
123 zm = x(3,n0) - z0
124
125
126 nc = nc + 1
127 iadll(nc+1) = iadll(nc) + 5
128 ik = iadll(nc)
129 lll(ik) = n0
130 jll(ik) = 1
131 sll(ik) = 0
132 xll(ik) =-one
133 ik = ik+1
134 lll(ik) = n1
135 jll(ik) = 1
136 sll(ik) = 0
137 xll(ik) = half
138 ik = ik+1
139 lll(ik) = n2
140 jll(ik) = 1
141 sll(ik) = 0
142 xll(ik) = half
143 ik = ik+1
144 lll(ik) = n0
145 jll(ik) = 5
146 sll(ik) = 0
147 xll(ik) = zm
148 ik = ik+1
149 lll(ik) = n0
150 jll(ik) = 6
151 sll(ik) = 0
152 xll(ik) =-ym
153
154
155 nc = nc + 1
156 iadll(nc+1) = iadll(nc) + 5
157 ik = iadll(nc)
158 lll(ik) = n0
159 jll(ik) = 2
160 sll(ik) = 0
161 xll(ik) =-one
162 ik = ik+1
163 lll(ik) = n1
164 jll(ik) = 2
165 sll(ik) = 0
166 xll(ik) = half
167 ik = ik+1
168 lll(ik) = n2
169 jll(ik) = 2
170 sll(ik) = 0
171 xll(ik) = half
172 ik = ik+1
173 lll(ik) = n0
174 jll(ik) = 6
175 sll(ik) = 0
176 xll(ik) = xm
177 ik = ik+1
178 lll(ik) = n0
179 jll(ik) = 4
180 sll(ik) = 0
181 xll(ik) =-zm
182
183
184 nc = nc + 1
185 iadll(nc+1) = iadll(nc) + 5
186 ik = iadll(nc)
187 lll(ik) = n0
188 jll(ik) = 3
189 sll(ik) = 0
190 xll(ik) =-one
191 ik = ik+1
192 lll(ik) = n1
193 jll(ik) = 3
194 sll(ik) = 0
195 xll(ik) = half
196 ik = ik+1
197 lll(ik) = n2
198 jll(ik) = 3
199 sll(ik) = 0
200 xll(ik) = half
201 ik = ik+1
202 lll(ik) = n0
203 jll(ik) = 4
204 sll(ik) = 0
205 xll(ik) = ym
206 ik = ik+1
207 lll(ik) = n0
208 jll(ik) = 5
209 sll(ik) = 0
210 xll(ik) =-xm
211
212
213 nc = nc + 1
214 iadll(nc+1) = iadll(nc) + 4
215 ik = iadll(nc)
216 lll(ik) = n1
217 jll(ik) = 1
218 sll(ik) = 0
219 xll(ik) =-half
220 ik = ik+1
221 lll(ik) = n2
222 jll(ik) = 1
223 sll(ik) = 0
224 xll(ik) = half
225 ik = ik+1
226 lll(ik) = n0
227 jll(ik) = 5
228 sll(ik) = 0
229 xll(ik) = zr
230 ik = ik+1
231 lll(ik) = n0
232 jll(ik) = 6
233 sll(ik) = 0
234 xll(ik) =-yr
235
236
237 nc = nc + 1
238 iadll(nc+1) = iadll(nc) + 4
239 ik = iadll(nc)
240 lll(ik) = n1
241 jll(ik) = 2
242 sll(ik) = 0
243 xll(ik) =-half
244 ik = ik+1
245 lll(ik) = n2
246 jll(ik) = 2
247 sll(ik) = 0
248 xll(ik) = 0.5
249 ik = ik+1
250 lll(ik) = n0
251 jll(ik) = 6
252 sll(ik) = 0
253 xll(ik) = xr
254 ik = ik+1
255 lll(ik) = n0
256 jll(ik) = 4
257 sll(ik) = 0
258 xll(ik) =-zr
259
260
261 nc = nc + 1
262 iadll(nc+1) = iadll(nc) + 4
263 ik = iadll(nc)
264 lll(ik) = n1
265 jll(ik) = 3
266 sll(ik) = 0
267 xll(ik) =-half
268 ik = ik+1
269 lll(ik) = n2
270 jll(ik) = 3
271 sll(ik) = 0
272 xll(ik) = half
273 ik = ik+1
274 lll(ik) = n0
275 jll(ik) = 4
276 sll(ik) = 0
277 xll(ik) = yr
278 ik = ik+1
279 lll(ik) = n0
280 jll(ik) = 5
281 sll(ik) = 0
282 xll(ik) =-xr
283
284
285
286 nc = nc + 1
287 iadll(nc+1) = iadll(nc) + 9
288 ik = iadll(nc)
289 lll(ik) = n1
290 jll(ik) = 4
291 sll(ik) = 0
292 xll(ik) =
alpha*xt1(1)
293 ik = ik+1
294 lll(ik) = n1
295 jll(ik) = 5
296 sll(ik) = 0
297 xll(ik) =
alpha*xt1(2)
298 ik = ik+1
299 lll(ik) = n1
300 jll(ik) = 6
301 sll(ik) = 0
302 xll(ik) =
alpha*xt1(3)
303
304 ik = ik+1
305 lll(ik) = n2
306 jll(ik) = 4
307 sll(ik) = 0
308 xll(ik) = xt2(1)
309 ik = ik+1
310 lll(ik) = n2
311 jll(ik) = 5
312 sll(ik) = 0
313 xll(ik) = xt2(2)
314 ik = ik+1
315 lll(ik) = n2
316 jll(ik) = 6
317 sll(ik) = 0
318 xll(ik) = xt2(3)
319
320 ik = ik+1
321 lll(ik) = n0
322 jll(ik) = 4
323 sll(ik) = 0
324 xll(ik) =-
alpha*xt1(1)-xt2(1)
325 ik = ik+1
326 lll(ik) = n0
327 jll(ik) = 5
328 sll(ik) = 0
329 xll(ik) =-
alpha*xt1(2)-xt2(2)
330 ik = ik+1
331 lll(ik) = n0
332 jll(ik) = 6
333 sll(ik) = 0
334 xll(ik) =-
alpha*xt1(3)-xt2(3)
335
336
337 nc = nc + 1
338 iadll(nc+1) = iadll(nc) + 6
339 ik = iadll(nc)
340 lll(ik) = n1
341 jll(ik) = 4
342 sll(ik) = 0
343 xll(ik) = yt1(1)
344 ik = ik+1
345 lll(ik) = n1
346 jll(ik) = 5
347 sll(ik) = 0
348 xll(ik) = yt1(2)
349 ik = ik+1
350 lll(ik) = n1
351 jll(ik) = 6
352 sll(ik) = 0
353 xll(ik) = yt1(3)
354 ik = ik+1
355 lll(ik) = n0
356 jll(ik) = 4
357 sll(ik) = 0
358 xll(ik) =-yt1(1)
359 ik = ik+1
360 lll(ik) = n0
361 jll(ik) = 5
362 sll(ik) = 0
363 xll(ik) =-yt1(2)
364 ik = ik+1
365 lll(ik) = n0
366 jll(ik) = 6
367 sll(ik) = 0
368 xll(ik) =-yt1(3)
369
370
371 nc = nc + 1
372 iadll(nc+1) = iadll(nc) + 6
373 ik = iadll(nc)
374 lll(ik) = n1
375 jll(ik) = 4
376 sll(ik) = 0
377 xll(ik) = zt1(1)
378 ik = ik+1
379 lll(ik) = n1
380 jll(ik) = 5
381 sll(ik) = 0
382 xll(ik) = zt1(2)
383 ik = ik+1
384 lll(ik) = n1
385 jll(ik) = 6
386 sll(ik) = 0
387 xll(ik) = zt1(3)
388 ik = ik+1
389 lll(ik) = n0
390 jll(ik) = 4
391 sll(ik) = 0
392 xll(ik) =-zt1(1)
393 ik = ik+1
394 lll(ik) = n0
395 jll(ik) = 5
396 sll(ik) = 0
397 xll(ik) =-zt1(2)
398 ik = ik+1
399 lll(ik) = n0
400 jll(ik) = 6
401 sll(ik) = 0
402 xll(ik) =-zt1(3)
403
404
405 nc = nc + 1
406 iadll(nc+1) = iadll(nc) + 6
407 ik = iadll(nc)
408 lll(ik) = n2
409 jll(ik) = 4
410 sll(ik) = 0
411 xll(ik) = yt2(1)
412 ik = ik+1
413 lll(ik) = n2
414 jll(ik) = 5
415 sll(ik) = 0
416 xll(ik) = yt2(2)
417 ik = ik+1
418 lll(ik) = n2
419 jll(ik) = 6
420 sll(ik) = 0
421 xll(ik) = yt2(3)
422 ik = ik+1
423 lll(ik) = n0
424 jll(ik) = 4
425 sll(ik) = 0
426 xll(ik) =-yt2(1)
427 ik = ik+1
428 lll(ik) = n0
429 jll(ik) = 5
430 sll(ik) = 0
431 xll(ik) =-yt2(2)
432 ik = ik+1
433 lll(ik) = n0
434 jll(ik) = 6
435 sll(ik) = 0
436 xll(ik) =-yt2(3)
437
438
439 nc = nc + 1
440 iadll(nc+1) = iadll(nc) + 6
441 ik = iadll(nc)
442 lll(ik) = n2
443 jll(ik) = 4
444 sll(ik) = 0
445 xll(ik) = zt2(1)
446 ik = ik+1
447 lll(ik) = n2
448 jll(ik) = 5
449 sll(ik) = 0
450 xll(ik) = zt2(2)
451 ik = ik+1
452 lll(ik) = n2
453 jll(ik) = 6
454 sll(ik) = 0
455 xll(ik) = zt2(3)
456 ik = ik+1
457 lll(ik) = n0
458 jll(ik) = 4
459 sll(ik) = 0
460 xll(ik) =-zt2(1)
461 ik = ik+1
462 lll(ik) = n0
463 jll(ik) = 5
464 sll(ik) = 0
465 xll(ik) =-zt2(2)
466 ik = ik+1
467 lll(ik) = n0
468 jll(ik) = 6
469 sll(ik) = 0
470 xll(ik) =-zt2(3)
471
472 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB